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University of Macau 
Abstract 

Numerical Investigation of Nonisothermal Reversed 
Stagnation-point Flow 
by Chio Chon Kit 
Thesis Supervisor: 
Associate Professor Sin Vai Kuong 
Electromechanical Engineering 

This thesis investigates the nature of the development of two-dimensional laminar 
nonisothermal flow of an incompressible fluid close to the reversed stagnation- 
point. Proudman and Johnson (1962) [Ij flrst studied the flow and obtained an 
asymptotic solution by neglecting the viscous terms. This is not practice in ne- 
glecting the viscous terms within the total flow fleld. Viscous terms in this analysis 
are now included, and two-dimensional nonisothermal reversed stagnation-point 
flow is investigated by solving the Navier-Stokes equations coupled to energy 
equation. 
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CHAPTER 1: INTRODUCTION 



The Navier-Stokes equations describe the motion of fluid substances by applying 
Newton's second law to fluid motion. It is wonder that given their wide range 
of practical uses, mathematicians are difficult or impossible to obtain an exact 
solution in almost every real situation because of the analytic difficulties associ- 
ated with the nonlinearity due to convective acceleration. The existence of exact 
solutions is fundamental not only in their own right as solutions of particular 
flows, but also are useful as accuracy checks for numerical solutions. 

Computational fluid dynamic modeling has been a very active area of research in 
recent years as evidenced by numerous papers in the literature. Advances in com- 
puter capacity concurrent with the maturation of flow and heat transfer modeling 
have made feasible these coupled simulations. The goal of research in this area is 
to make simulations simple in the design and analysis environment for real- world 
applications. Such capability would be very beneflcial to those industries, includ- 
ing enhanced oil recovery, which is a technique for increasing the amount of crude 
oil that can be extracted from an oil fleld. 

1.1 PREVIOUS AND RELATED WORK 

Several application of such have also appeared in the recent literature, for exam- 
ple, in some simplifled cases a fluid travels through a rigid body (e.g., missile. 
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sports ball, automobile, spaceflight vehicle), or in oil recovery industry crude oil 
that can be extracted from an oil fleld is achieved by hot water injection, as 



shown in Fig. (1.1), or equivalently, an external flow impinges on a stationary 
point called stagnation-point that is on the surface of a submerged body in a flow, 
of which the velocity at the surface of the submerged object is zero. Moreover, 
the streamline is perpendicular to the surface of the rigid body. The study of the 
flow motion at stagnation point is of importance in oil recovery industry that de- 
velops techniques to efliciently recover oil, gas, and other minerals while reducing 
environmental impacts using various pollution remediation and greenhouse gas 
reduction techniques. 



CO2 
Injector well 

I 



Producer well 
> 




Figure 1.1: Oil recovery industry 



The classic problems of two-dimensional stagnation-point flows can be analyzed 
exactly by Hiemenz [5j. The result is an exact solution for flow directed perpen- 
dicular to an inflnite flat plate. Howarth [6j and Davey [7J extended the two- 
dimensional and axisymmetric flows to three dimensions, and Wang [8j studied 
the case for obliquely-impacting jets. The similarity solutions for the tempera- 
ture fleld were studied by Eckert [9j. Case corresponding a step change in wall 
temperature or in wall heat flux in laminar steady flows at a stagnation point 
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has been also investigated by several authors (see Chao et al. [10], Sano [H] 
and Gorla [12j). Further, Lok et al. [13j investigated the mixed convection near 
non-orthogonal stagnation point flow on a vertical plate with uniform surface 
heat flux, where the results published are very good with present value of the 
normalized temperature at the wall for the constant wall temperature boundary 
condition. 



On the contrary, when the external flow is extracted away from the stagnation- 



point shown in Fig. (1.1), the flow in the vicinity of this "reversed stagnation 
point" is governed by boundary-layer separation and vorticity generation and the 
reversed stagnation-point flow develops. Reversed stagnation-point flow is a flow 
in which the component of velocity normal to a wall is outward the wall every- 
where in the region concerned. Reversed stagnation-point flows against an inflnite 
flat wall do not have analytic solution in two dimensions, but certain reverse flows 
have solution in three dimensions [7j. 



Proudman and Johnson [T] first suggested that the convection terms dominate in 
considering the inviscid equation in the body of the fluid. By introducing a very 
simple function of a particular similarity variable and neglecting the viscous forces 
in their analytic result for region suflicient far from the wall, they obtained an 
asymptotic solution in reversed stagnation-point flow, describing the development 
of the region of separated flow for large time t. In their solution, the phenomenon 
of separation is described near a plane that represents the rear-stagnation point 
of a cylinder is set in motion impulsively with a constant velocity normal to 
the surface of the plane. Robins and Howarth [14] have recently extended the 
asymptotic solution, flnding higher order terms by singular perturbation methods. 
They indicated that the viscous forces cannot be ignored in the governing equation 
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because of a consistent asymptotic expansion in both this outer inviscid region 
and also in the inner region near the plane. Smith [15j generalized the solution 
of Proudman and Johnson with both viscous and convection terms in balance by 
considering the monotonic potential flow when the time is relatively large. Shapiro 
[16] obtained a solution for unsteady reversed stagnation-point flow with injection 
or suction. These unsteady flows flt within a class of similarity transformations 
originally identifled by Birkhoff using a group-theoretic approach. 

1.2 THESIS OBJECTIVE 

This thesis focuses on the challenging problem of numerical modeling for a non- 
isothermal reversed stagnation-point flow. The primary objective of the present 
study is to determine the main characteristics of the flow at reversed stagnation- 
point. This includes the flow proflle, the separation zone, the dividing streamline 
and the nonisothermal temperature proflle. The other objective is to briefly verify 
the simulation data for the reversed stagnation point. 

1.3 THESIS OUTLINE 

Chapter 2 is dedicated to developing the governing equation of flow and heat 
transfer modeling used in this research. Chapter 3 investigates the nature of 
the development of two-dimensional laminar flow of an incompressible fluid near 
the reversed stagnation-point. Similarity solutions of two-dimensional reversed 
stagnation-point flow are investigated by simplifying the full Navier-Stokes equa- 
tions coupled to the energy equation, describing the motion of nonisothermal fluid 
substances. The model is valid if the fluid velocity is small compared with the 
speed of sound and the fluid is treated as Newtonian. Chapter 4 describes the 
implementation of the algorithm used to achieve the objectives outlined above. 
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Specific information is provided about the individual solvers and details of the in- 
terface of simulation. Details of the mesh generation and the efficacy of the CFD 
solver are provided. Chapter 5 focuses on numerical simulations on the reversed 
stagnation-point flow. This chapter provides result and discussion of the model 
provided. The final chapter summarizes and provides conclusions of the research, 
and recommends future work. 



CHAPTER 2: GOVERNING EQUATIONS 



The nonlinear behavior of fluid flow is emphasized in this chapter. The Navier- 
Stokes equations describe the motion of a fluid in two- or three-dimensional space. 
These equations are to be solved for an unknown velocity vector and pressure. We 
restrict attention here to the reversed stagnation-point flow in an incompressible 
fluids domain. To construct an effective method for handling the Navier-Stokes 
equations, which are systems of partial differential equations, a similarity trans- 
formation is applied and a simplifled similarity equation is considered. 

2.1 INCOMPRESSIBILITY 

In an incompressible fluid, the density of an element of fluid is not affected by any 
changes in pressure. If the relative speeds within a flow are low enough (typically 
Mach number less than 0.3), thermodynamic effects and density changes due to 
changes in pressure become negligible. If density is constant and mass is conserved 
so is volume. This condition is called the equation of continuity and expressed 
mathematically as the divergence of flow velocity V is zero 

V-V = (2.1) 

Essentially what goes into a differential volume must exit it simultaneously. Cou- 
pling this equation with conservation of momentum makes the system fully de- 
termined, without need of the energy equation or an equation of state, and yields 
extremely efficient simulations. The flows and solution methods can be greatly 
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different; however, they aU start with the same underlying defined as the diflFer- 
ential element of the continuous Navier-Stokes equations. 

2.2 MOMENTUM EQUATIONS 

Models for Newtonian fiuids undergoing incompressible fiow make use of the ap- 
proximation that dynamic viscosity /x is a constant. Performing a force balance 
and making use of the continuity equation leads to the Navier-Stokes equations 

p^ = -Vp + ^N^V + f (2.2) 

DV / Dt is the material derivative of flow velocity 

DV dV ^ ^ 

_ = _ + (V.V)V (2.3) 

representing the convective acceleration in the fiuid motion. The physical princi- 



ple of momentum transfer is Newton's second law. Equation (2.2) is just Newton's 
second law F = ma for a fiuid element subject to the external force / and to the 
forces arising from pressure gradient — and viscosity /xV^V. 

2.3 ANALYTICAL ANALYSIS 

We begin with writing the governing equation in conservative velocity form in the 

Cartesian coordinates [18j and neglecting the external force /: 

du dv ^ , 

dx dy 

du du du 1 dp f d^u d^u\ 

dt dx dy pdx \dx^ dy^ J 



dv ^ dv ^ dv 1 dp ^ / d'^v ^ d'^v . (2 4) 

dt dx dy pdy dy"^ J 
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Here u and v are the components of flow velocity V{u^v)^ p is the fluid density, 
p is the fluid pressure, u = p/p is the kinematic viscosity. The viscous fluid flows 



in a rectangular Cartesian coordinates (x, z). Fig. 2.1, illustrates the motion of 
external flow directly moving perpendicular out of an inflnite flat plane wall. The 
origin is the so-called stagnation point and z is the normal to the plane. 



> 


V 






u = 


Tu, 

i > 



X 



Figure 2.1: Coordinate system of nonisothermal reversed stagnation-point flow 

What we are concerned about is the two-dimensional reversed stagnation-point 
flow in unsteady state. The total fluid domain is bounded by an inflnite plane 
^ = 0, the fluid remains at rest when time t < 0. At t = 0, it starts impulsively 
in motion which is determined by the stream function 



= —axy 



(2.5) 



where a is a positive rate of strain [T8]. At large distances far above the planar 
boundary, the existence of the potential flow implies an inviscid boundary condi- 
tion. Far away from the wall the flow is of a constant Vq along the y-axis. Because 
of the axisymmetric conflguration, the flow fleld is considered in right-hand side 
region only. For such a flow the components of velocity are easily given from the 
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relationships 

u —ax (2.6a) 
v = Vo (2.6b) 

Here a is a constant proportional to Vq/L^ Vq is the external flow velocity remov- 
ing from the plane and L is the characteristic length. We have 2i = at x = 
and ^ = at = 0, but the no-slip boundary at wall {y = 0) cannot be satisfled. 



The equation of continuity (2.4a) is integrated by introducing the stream function 

u = — and V = — — (2.7) 
oy ox 

For reversed stagnation flow without friction (ideal fluid flow) , the stream function 
may be written as 

^ = ^l^^ = -A,xy (2.8) 
where is a constant and from which 

Ui = —AiX and Vi = A^y. (2.9) 

We have = at x = and = at ?/ = 0, but the no-slip boundary at wall 
{y = 0) cannot be satisfled. 

Since for a (real) viscous fluid the flow motion is determined by only two factors, 
the kinematic viscosity u and a, consistent with the initial and boundary condi- 
tions that ip is proportional to x for all value of y and t. Provided that the surface 
is an inflnite plane wall, a following modifled stream function is introduced, see 
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Proudman and Johnson 



^ = -^/Apxf{r],T) (2.10a) 

[a 

r]=J-y (2.10b) 

T = At (2.10c) 



where i] is the non-dimensional distance from waU and r is the non-dimensional 
time. Noting that the stream function automatically satisfies equation of continu- 



ity (2.4a). Substituting u and v into the governing equations results a simplified 



partial differential equation. From the definition of the stream function, we have 

« = ^ = -Axf^ (2.11a) 
v = -^ = ^f (2.11b) 

Note that A has the dimension as "1/time". The governing equations can be 
simphfied by a similarity transformation when several independent variables ap- 
pear in specific combinations, in flow geometries involving inflnite or semi-inflnite 
surfaces. By introducing coordinate variable transformation, the number of in- 
dependent variables is reduced by one or more. The original system of partial 
differential equations can be simplifled into the following pair of partial differen- 
tial equations 



—A xfrjT + A x{frj) — A xffrjrj = ~ ^ ^frjijr] (2.12a) 



AVI^fr + A^//^ = + ^^/w (2.12b) 

The pressure gradient can be again reduced by a further differentiation equa- 



tion (2.12b) with respect to x. That is 
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and equation (2.12a) reduces to 

[fvr - Uv? + ffm - fvvvh = 0- (2.14) 
or the equation becomes a differential equation for / [19] 

fvT-ifv) + //r/r? - /7?7?7? = function of T only. (2.15) 
The initial and boundary conditions are 

f{v,0)=v iv^O) (2.16a) 

/(0,r) = /^(0,r) = (t^O) (2.16b) 
/(oo,r)-?7 (2.16c) 

These foUow from the impermeabihty condition of the waU (from 'u(x,0,r) = 
it foUows that /(0,r) = 0) and from the no-shp condition (from 2i(x,0,r) = it 
foUows that /ry(0,r) = 0). 

Proudman and Johnson suggested that at large distances from the wall {rj ^ oo) 
the velocity v^x^rj^r) should pass over smoothly into that for inviscid Vq. Here 
they have employed frjr]{oo) = frjr]r]{^) = 0, which implies that the flow matches 
smoothly with the inviscid flow as ry ^ (X). This leads to the condition /(oc, r) ^ 
7y, and thus, the last condition reduces the diflFerential equation ( 2.14D for / [19j 



frjr iff]) + ffrjT] frjrjr] — 1^ (2-17) 
with the boundary conditions 

/(0,r) = /^(0,r) = (2.18a) 

/^(oo,r) = l. (2.18b) 

Here frjrjr] is proportional to the viscous stress, frjr] is proportional to the shear 
stress, fri is proportional to the x-component of velocity in boundary layer and / 



is s proportional to the stream function. 



It should be noted that the dimensionless velocity distribution ff, is, from (2.11), 



independent of the length x, and thus equation (2.17) is a similarity equation 
of the full Navier-Stokes equations at two-dimension reversed stagnation-point. 
The coordinates x and y are replaced by a dimensionless variable rj. Under the 
boundary conditions /7^(oo,r) = 1, when the flow is in steady state such that 
f^j- = 0, the diflFerential equation has no solution. 



2.4 ENERGY TRANSPORT 

In this section our considerations of reversed stagnation-point flow until now have 
referred only to velocity fleld. Now we shall extend to include the temperature 
field in the nonisothermal flow which is at a temperature T different from that of 
the wall T^. It will be assumed that heat energy is transferred to the flow through 
the wall. Once the velocities are known from the flow analysis, the temperature 
distributions can be determined by solving the energy equation in the reversed 
stagnation-point flow. 



To include the temperature T in our analysis we must now turn to the thermo- 
dynamic properties of fluids. The principle of conservation of energy yields the 
thermal energy equation [20j for constant-property fluid: 



dT dT dT\ , fd^T d^T 

\- v. \- V I = I ^ A ^ 



dt ~^ ^ dx ~^ ^ dy 



+ 



\ dy 



(2.19) 



where k is the thermal conductivity, Cp is the heat capacity, and $ is deflned as 
2 /n..\2l 



dx J \dy 



dy dx J ?> \dx dy J 



(2.20) 
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and is called the viscous dissipation since it represents the irreversible conserva- 
tion of mechanical forms of energy to a thermal form. 

If both velocity field and temperature field exist, there is generally also a coupling 
between these two fields. Since the velocity components u and v appear in the 
energy equation, a simplification of the energy equation requires to know the 
actual value of the velocity components. This velocity field would be identical to 
the velocity components in the reversed stagnation-point fiow 

u = -Axfrj (2.21a) 
V = VAi^f (2.21b) 



To transform equation (2.19) into a nondimensional form, it is convenient to work 



with a dimensionless temperature 9 [18j: 

0{v,t) = ^~^;: (2.22) 

where and Too are the wall temperature and ambient temperature. Consider- 
ing the case that both and Too are constant, the required boundary conditions 
are 

T(0,t) = T^, T(oc,t) =Too. (2.23) 

The fiuid temperature T can be treated as a function of rj and r only. Under the 
assumption that the viscous dissipation is negligible compared to conduction at 
the wall, we may write the energy equation in the form 



(2.24) 



subject to the boundary conditions 

e(0,r) = ^((X),r) = l (2.25) 
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Equation ( |2.24D is a second-order partial differential equation with variable coef- 
ficient f{r]^r) and the Prandtl number Pr = pcpu/k is assumed to be constant. 
Consider the fluid of which Pr = 1, the thermal boundary layer and the velocity 



boundary layer collapse, and thus, substituting 6 = f \ equation (2.17) and (2.24) 
represent the same equation. It is noticed that in these nonisothermal flows the 
velocity fleld is decoupled from the temperature fleld if the kinematic viscosity v 
is constant and is assumed to be independent of the temperature and pressure. 
This assumption is valid when the temperature and pressure differences are small 
within the boundary layer. 



CHAPTER 3: FLOW ANALYSIS 



We complete the governing equations of viscous reversed stagnation-point flow 
by discussing similar flow. Our objective is to obtain a similarity solution of 
the governing equation. Generally speaking, a similarity solution is one in which 
the number of variables can be reduced by a coordinate transformation. Let us 
discuss various laminar similarity solutions. 



3.1 INVISCID SOLUTION 

Proudman and Johnson P first thought over the early stages of the diffusion 
of the initial vortex sheet at ^ = 0. The idea was to divide the flow into two 
regions: an outer flow region that is inviscid and can sometimes be approximated 
as potential flow, and an inner flow region where the viscous forces are of the 
same order as the inert ial forces. The general feature of the predicted streamline 



pattern is sketched in Fig. (3.1). 



Proudman and Johnson suggested that, when the flow is near the wall region, the 
viscous forces are dominant, and the viscous term in the governing Navier-Stokes 
equations is important only near the boundary. On the contrary, the viscous 
forces were neglected far away from the wall. The convection terms dominate the 
motion of external flow in considering the inviscid equation in the fluid. Ignoring 
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X 



Figure 3.1: Streamlines of reversed stagnation-point flow 



the viscous stress in equation (2.17) yields an inviscid equation 



fvr-{fvf + ffvV + '^ = 0- 



(3.1) 



They contemplated the similarity of the inviscid equation in the form 



/(77,r) = A(r)F(7), 7 = VA(r) 



(3.2) 



Substituting equation (3.2) in (3.1) results in 

A 



A 



7F" 



F'2 + pp" 



■1, 



SO that 



Y constant = fc, or A = e^^ 

A 



(3.3) 



(3.4) 



A solution to this equation that satisfies F = 1 with exponential error as 77 ^ cxo is 
only possible when k = 1; Proudman and Johnson finally obtained an asymptotic 
similarity solution of 

2e^ 



(3.5) 



where c is a constant of integration, representing the uncertainty in the precise 



position of the time origin. In the asymptotic solution (3.5), the constant c always 
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appears multiplying the similarity variable, i.e. cr]e~^ . A change from r to r + Ar 
can be included in the constant c. The improved numerical evaluations of Robins 
and Howarth [14j estimated the value of c to be 3.51. This solution describes an 
exponential decay of vorticity in the outer region, moving away from the plane 
with a constant velocity. 



3.2 VISCOUS SOLUTION 

The viscous layer develops as a consequence of the no-slip boundary condition at 
the wall. In the inner region the viscous term cannot be neglected and a further 
solution must be found which satisfies the no-slip condition on the wall. When 



r ^ oc, the solution (3.5) yields the steady fiow 

/ - -77 and / - -1 (3.6) 
which becomes the outer boundary condition for the viscous fiow near the bound- 



ary. Substituting in equation (2.17) yields 

f" - ff" + if? -1 = (3.7a) 
/(O) = /(O) = (3.7b) 
/(oo) = -1 (3.7c) 

This is exactly the classic stagnation-point problem (Hiemenz [5j) by changing 
the sign in /. It is a third-order nonlinear ordinary differential equation and does 
not have an analytic solution, and thus it is necessary to solve it numerically. The 



numerical solution of classic stagnation-point problem is shown in Figure (3.2). 



Although an asymptotic solution was obtained, it can easily been observed that 
this is not valid when the viscous term frjrjrj is neglected within the total fiow field. 
Robins and Howarth [14j studied higher order terms by singular perturbation 
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Figure 3.2: Numerical solution of classic stagnation-point problem 



methods and indicated that a consistent asymptotic expansion occurs in both 
outer inviscid region and also in the inner region that must exist close to the wall 
where the viscous forces need to be included. It is not quite appropriate to say 
Proudman and Johnson are wrong because of neglecting the viscous term in their 
analytic result for region sufficient far from the wall, but in such a case neglecting 
the viscous terms within the total flow fleld can be improved. The next section 
will discuss the nonexistence of the exact solutions with the boundary condition 
/7^(oc,r) = 1 for steady case. 



3.3 INSOLUBILITY IN STEADY STATE 

In this section, a mathematical proof indicates that all of the steady solutions, 
however, do not satisfy the boundary condition fj^^oo^r) = 1. 
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When fnr^O, Eq. (|2.17p reduces to 



/" - //" + (/)' -1 = 

/(O) = /(O) = 
/(oo) = 1 

where the prime denotes the derivative with respect to rj. 



(3.8a) 
(3.8b) 
(3.8c) 



Lemma 1 No solution f'ii]) exists which has stationary value of 1 for finite rj. 



Proof. Rearrange equation (3.8) yields 



= 1 - ifr + ff 



f\2 , r nil 



(3.9) 



Suppose for rj = tjq^ we have f^rjo) = 1 and f^^rjo) = 0. Afterwards, it follows 
from the derivatives of equation (3.9) that /^^^ and all higher derivatives are zero 



when 7] = r]Q. Consider a variable transformation 

zu{r]o) = 1 

Expand the function into Taylor's series near rjQ, we have 



(3.10) 



OO 



n=0 



oo 



n=l 



nl 



n=l 



nl 



= 1 



Hence, the boundary condition f\0) = is thus not satisfied and the Lemma is 
proved. 
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Lemma 2 When has a stationary value, if \f^\ < 1 it is a minimum and if 
\f^\ > 1 it is a maximum. 



Proof: From equation (3.9), when has a stationary value, it means /^^ = and 
equation (3.9) becomes 

/" = 1 - (/)' (3.11) 

If l/^l < 1, /^^^ > and it is minima. Else if |/^| > 1, /^^^ < and it is maxima. 
Eventually, the lemma is proved. 

Lemma 3 // f^\r]) vanishes for rj = 771,772 ^'^th 771 < 772 < then the 
sequence f'{r]i) does not tend to 1 as rji ^ 00. 

Proof: Consider a region where (771,772) is far away from the origin. Multiply /^^ 



to equation (3.9) and integrate it between 771 and 772 with respect to 77. 

f"f" = f" - if'ff + fiff 
'f"f"dv= / [f"-{f?f"+f{f"f]d^ 



2 

^ ^ Jr]i 

When f\oo) 1, it is required that f^\r]i) = f^\r]2) = and thus 

[f'-l{f')X = -L 

whereas L = f {f ) drj is always positive and we can obtain 
Jr]l 



[/-^(/?]??<o 



f'{m)-l[f'{v2)?<f'{m)-l[f'{vi)] 

[f'{V2)f - S/fe) > [f'{vi)f - 3/(r7i) 
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Consider G = f'^ - 3/' function of then 



= 3/'^ - 3 



As = 1, then G\l) = 0, which makes G a minimum. We do not have f\rij) = 1 
as ry^ ^ oo 



Theorem 1 Given any f\r]) 1 as rj ^ oo, no solution of equation (3.8) exists. 



Proof : When \ f'\ < 1, since ^ 1 as ^ oc, then /^^ must be greater than 



zero. Hence, recaU from equation (3.9), 



/" = 1 - iff + ff > 0. 
for aU T] > r]Q. After integrating f"'{r]) > from tjq to ij > tjq^ we have 

f'iv) > fivo) = K>0. 
Another integration from tjq to r] > tjq yields 

f'iv)>f'ivo) + Kiv-m)- 

By Lemma 2, f\r]) has at most one stationary value because one cannot have 
two consecutive stationary values which are both minima. Since f^\r]) > 0, when 
7y ^ oc, f\r]) oo. It violates that f\r]) ^ 1. A similar argument shows that a 
solution cannot approach to 1 when |/^| > 1. 



The remaining option is that f{r]) oscillates about 1 as f\oo) 1. But this 
would imply an infinite number of changes in concavity of as 77 ^ oc. Once 
f\r]) vibrates from concave upward (/^^^ > 0) to concave downward (/^^^ < 0), as 
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f{r]) becomes great enough, it cannot turn concave upward again [21]. It requires 
a point rji such that 

/(4)(^l) = 



and 



Differentiating equation (3.9) gives 



(3.12) 



and 



/(5) = ffii) _ ^ff 



Evaluating these at rji yields 



(3.13) 



If f"{r]i) 7^ we get an immediate contradiction. On the other hand, if f"{r]i) 



0, equation (3.12) implies that /^^^ = 0. But from equation (3.9) we see that 
f^\vi) = f''\vi) = implies the desired contradiction that 

f{r,) = 1. 

Thus, since f\r]) cannot ultimately approach to 1 from above or below, nor in 



an oscillatory manner, no solution to equation (3.8) exists. Reversed stagnation- 



point flow against an impermeable flat wall does not exist in two-dimensional 
steady case. 



3.4 FINITE-DIFFERENCE FORMULATIONS 



Similarity solutions of reversed stagnation-point flow with different boundary con- 
ditions have been published in [22l [23l [24]. Numerical simulation of reversed 
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stagnation-point flow with full Navier-Stokes equations has been studied in [25] . 
According to the previous work, the governing equations in reversed stagnation- 
point flow are 



Or^ri - Pr fOr^ = Pr Or (2.24) 



The above equations subject to the boundary conditions (2.17) and (2.24) are 
nonlinear third-order partial differential equations. They do not admit similarity 
solution and numerical or perturbation methods are required to solve the problem. 

We shall, however, use here a numerical method. It is an implicit flnite-diflFerence 
method with second-order accuracy. The partial diflFerential equations can be ex- 
pressed as approximate expressions, so that it is easy to program the solution of 
large numbers of coupled equation. 

We start with rewriting the partial differential equations in the form: 



fr]r — fr]r]r] + (///) 1 + f fijr] (3.14a) 
1 

Pr 



^r = ^Jm-fOr^ (3.14b) 



and introducing the new dependent variables 



h = l- (3.15a) 
g = e (3.15b) 



The equations can be rewritten as 

hr = hnr, + 2h-h^ + hr, j {l-h) dr) (3.16a) 
9t = ^9r]r] - 97] J {I - h) drj (3.16b) 



25 



We now contemplate the net rectangle in the t — rj plane shown in Fig. (3.3) and 
the net points defined as below: 

rjO = 0, rij = rij_i + Ar], j = 1, 2, ...J, rjj = ?]oo 



r 



= 0, r" = T^-^ + +Ar, n=l,2,...J, 



Here n and j are just the sequence of numbers that indicate the coordinate loca- 
tion, not tensor indices or exponents. The partial differential equations are easily 

V 



Vj+i 

Vj - 
Vj-i -- 



Pi 



Arj 



O 



At 
H h 



Figure 3.3: Net rectangle for finite-difference method 

discretized by central difference representations with second-order accuracy, for 
example the finite diflFerence forms for any points are 



hr 



and 



hr 



2Ar] 



K^l - + K-1 



(3.17) 



(3.18) 



When i = 0, since the value of hf_i is not logical, the derivative is replaced by 
the forward diflFerence with second-order accuracy 
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The finite-difference form of the ODE is written at the midpoint (r^, rjj)^ the 
discretized equation takes the form 

-K_ - 2hr' + hti 



At 



(K? - (1 - h) dv (3.20a) 



9i 9i ^ 9i^i ^ 9i-i 9i^i 9i-i 

At 



iAr] 



Pr{Arif 2Ar) Jq 

This procedure yields the foUowing hnear tridiagonal system: 



(1 - h) dr] (3.20b) 



n+1 Rh^^^ 



un _-Ln I 



Ed - 





9? - 



(3.21a) 



5^(1 - hf) (3.21b) 




where /3 = Ar/(A7y)^. 

The initial conditions are the solutions of the following second-order linear parabolic 
differential equations 



h 
9r 



1 



Pr 



9r]r] 



(3.22a) 
(3.22b) 



As can be seen from the energy equation (2.24), equations (3.22) are identical 



to the heat conduction equation for one-dimensional unsteady temperature field, 
and thus, there are many solutions to these differential equations in [26]. The 
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desired solutions of (3.22) have the form 



/i = 1 - erf i^-^j (3.23a) 

^ = erf ( — I (3.23b) 

where the error function erf(z) is defined as 

erf(z) = ^ rexp{-C'^)di (3.24) 
V^r Jo 

When r ^ 0, the boundary conditions are convenient to write in the form 

Ag = sS=0, '.» = erf(^), = erf (^-^j (3.25) 



Equations (3.21) are defined as being imphcit, as more than one unknown appears 
in the left hand side. They are unconditionaUy stable, however, set of linear al- 
gebraic equations is required to be solved by the tridiagonal matrix algorithm 
(TDM A), also known as the Thomas algorithm, which is a simplified form of 
Gaussian elimination that is applied to evaluate tridiagonal systems of equations. 

For the stability of the diffusion difference equation, the condition of /5 < ^ must 
be satisfied. The procedure is straightforward, except for the algebra. The result- 
ing algorithm of the finite-difference method is written in MATLAB, a numerical 
computing environment allowing matrix manipulations and plotting of functions 
and data. At our level of discretization, however, we are only able to resolve in 



small time range. The numerical results of (3.21) are presented in Figures (3.4) 
to (ra. 



From the Proudman- Johnson solution (3.5), we have 

log h = —crje^ + log 2, 



(3.26) 
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so the graph of log/i against rje^ should provide a straight hne of gradient —c 
if the Proudman- Johnson solution holds. In Figure ( |3.4| 3), the graph of log/i 
against e~^^^'^ is plotted for different values of r. As can be seen, the parallel 
straight lines for large values of r agrees well for the Proudman- Johnson solution. 
The value of c calculated from the gradient of the straight line in Figure ( |3.4| 3) is 
c = 3.5, which agrees to the estimation of Robins and Howarth. 



The following pages (Figures. (3.5) to ( |3.6[ )) show the numerical solution of tem- 
perature distributions with Pr. It is noted that the dimensionless wall tempera- 
ture gradient g\0) raises with increase of Prandtl number, but the temperature 
boundary layer thickness decrease with increase of Prandtl number. Prandtl num- 
ber is the characteristic number for thermal boundary layers and heat transfer in 
forced convection. It can be explained by the definition of Prandtl number that 
inversely proportional to the thermal diffusivity a. Prandtl number is a ratio of 
two quantities which characterize the momentum and heat transport of fiuid. In 
heat transfer problems, the Prandtl number controls the relative thickness of the 
momentum and thermal boundary layers. When Pr is small, it means that the 
heat diffuses very quickly compared to the velocity field. This means that for 
liquid metals the thickness of the thermal boundary layer is much bigger than the 
velocity boundary layer. 



Figure (3.7) shows the pressure distribution along the x- and ry-direction at dif- 
ferent values of r = 3,4 and 5. Lines without markers denote results obtained 
from numerical procedure and dotted lines are from asymptotic solution. Far 
away from the wall region, the solution agrees remarkably well for smaller values 
of r with the known asymptotic solution, thus confirming the predictions of the 
analytical solution. On the other hand, discrepancy occurs as a larger value of 




Figure 3.4: Numerical Solution of equation (3.5) against (a) r, (b) ije 




(b) Pr = 1 

Figure 3.5: Asymptotic temperature solution g for various value of 7 




(b) Pr = 10 

Figure 3.6: Asymptotic temperature solution g for various value of 7 



32 



is applied in the numerical simulation. However, in the region near the stagna- 
tion point, a large difference is observed from the results obtained by these two 
methods in the region near the stagnation point. Numerical findings show that 
pressure profiles obtained from asymptotic solution and numerical simulation are 
in tremendously good agreement for smaller value of r. Discrepancy of results in 
pressure profiles increases for larger value of r. 




(c) r = 5 

Figure 3.7: Numerical relations of pressure profiles 
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CHAPTER 4: NUMERICAL SIMULATION 



As noted earlier, the simulation data at the reversed stagnation point is studied 
by solving the full Navier-Stokes and energy equations numerically, with the aid of 
a free, open source CFD software package of OpenFOAM. This chapter describes 
the fundamentals of the finite volume discretization. The technique has been 
described by many authors [171 EH EJ EHJ [29] and is applied to solve the reversed 
stagnation-point fiow. 

4.1 FINITE VOLUME METHOD 

The finite volume method (FVM) is a numerical technique that evaluates partial 
differential equations (PDEs) in the form of algebraic equations. Similar to the 
finite difference method, values are calculated at discrete places on a meshed ge- 
ometry. An advantage of the finite volume method is that it is easy to develop to 
enable unstructured meshes. This feature gives convenience in processing com- 
plicated geometries. Once the mesh of the domain is formulated, those governing 
equations are able to be solved. The method is applied in many computational 
fiuid dynamics packages, for instance, STAR-CD, FLUENT and OpenFOAM. 



When solving the Navier-Stokes equations, discretization of the solution domain 



is shown in Figure (4.1). The solution domain consists of a space and a time 



domain. The space domain is subdivided into a set of very small but finite-sized 
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cells or volumes covering the whole domain. Each cell is stored a set of governing 
equations to describe the physical phenomenon. Discretization of time involves 
in subdividing the domain into a set of time steps At, which may alter during a 
numerical simulation depending on some condition calculated during the simula- 
tion. 




At 

Time domain 



Figure 4.1: Discretization of the solution domain from [2] 
Discretization of space involves in subdividing the domain into a number of cells, 



or control volumes. A typical cell is shown in Figure (4.2). As can be seen the 
cells fill the computational domain without overlap. Dependent variables and 
other properties are principally assigned at the cell centroid, but it is possible to 
assign them on faces or vertices. The cell is bounded by a set of fiat faces, given 
the generic label. 
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Figure 4.2: Parameters in finite volume discretization from [2j 

4.2 FUNDAMENTAL EQUATIONS 

The purpose of equation discretization is to transform one or more governing 
equations into a corresponding system of algebraic equations. The solution of 
this system approximates the solution to the original partial differential equa- 
tions at certain locations in space and time. As for compressible flows, the mass 
conservation is a transport equation for density. With an additional energy equa- 
tion p can be constructed from a thermodynamic relation (ideal gas law) . 



And for incompressible flows, density variation is not correlated to the pressure 
field. Mass conservation is a constraint on the velocity field. Combined with the 
momentum equation, an equation for the pressure can be derived analytically. 
Consider the continuity and Navier-Stokes equations 

V-V = (4.1) 

+ V 'VV = --Vp + vV^V (4.2) 
ot p 
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where p is the density, p is the pressure, u is the kinematic viscosity. The principle 
of conservation of energy wiU yield the equation of energy for negligible viscous 
dissipation: 

— + V • ry = av^r (4.3) 

ot 

where a = k/pcp is thermal diffusivity. 



4.2.1 INTEGRAL FORM OF THE EQUATION 

If y is a closed region in space enclosed by a surface S, then 



V-TdV= / n-TdS 
V Js 

where n is the outward normal surface vector and F can represent any tensor 
field. The finite volume method performs well on the physical conservation and 



is adopted in the present study. Taking the volume integral on equations (4.1) 



and (4.2) and transforming equations into the surface integral forms using Gauss 



divergence theorem, we can get the integral form of the equations as follows 

n . VdS = (4.4) 

s 



^ I VdV+ f n-VVdS = -- [ n ■ pdS + u f n-VVdS (4.5) 
dtJv Js pJs Js 

4- f TdV+ f n-TVdS = a f n-VTdS (4.6) 
ot Jv Js Js 

4.2.2 FINITE VOLUME DISCRETIZATION 

When solving the Navier-Stokes equations, the region is often discretized using a 
staggered grid, in which the different unknown variables are not located at the 
same grid points. In the grid we shall use, the pressure p is located in the cell 
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centers, the horizontal velocity u in the midpoints of the vertical cell edges, and 
the vertical velocity v in the midpoints of the horizontal cell edges. 
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Figure 4.3: Staggered grid with boundary cells 

The computational domain is divided into a set of discrete volumes which do 
not overlap and fill the computational domain completely. The above equations 
are then volume-integrated over each individual finite volume. To convert the 
divergence terms into surface-integrated fiux terms. Gauss's theorem is used to 
reduce the problem. The divergence term is discretized to one of finding difference 
approximations for the fiuxes at the surface of the control volume based on the 
known cell-center values. The temporal derivatives can be discretized using finite- 
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difference approximations. The integrals can be replaced in the sum terms: 



fidS = s 



f 



I ndS = ^S 



n 



VdS = Y^S -V = ^F 



V 



f f 

dV = Vp 



(4.7a) 
(4.7b) 

(4.7c) 

(4.7d) 



where / is one face in the polyhedral cells and P represent the cell. Equations 



(4.4) and (4.5) are descretized as follows, 



f 



(4. 



(4.9) 



d_ 
dt 



TpVp + ^FT^ = a ^{S ■ V)TdS 



(4.10) 



/ 



/ 



Equations (4.8) to (4.10) are linearized by fixing the flux F because of the linear 



of the variable. As a consequence the vector equation can be decomposed into 



three component equations. From the equation (4.9), the discretized V can be 
expressed: 



Vp 



H{V) 

ap 



ap 



-Vp 



and also 



my) 

ap 



—Vp 

ap yf 



(4.11) 



(4.12) 



It is noticed that equation (4.9) divided by the finite volume is then converted to 
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the equation (4.11) but the pressure field is not discretized. The N-S equation is 
dependent on the pressure through the pressure gradient term in the momentum 
equation but we do not have a dependent pressure equation. If the fiow is com- 
pressible the continuity equation can be used to obtain the density field which 
can be applied to solve the pressure from an equation of state. 

On the contrary, for incompressible fiows, the continuity equation becomes an 
additional constraint on the velocity field. One way to overcome this difficulty is 
to build up a pressure field such that velocity satisfies the continuity equation. 



From the continuity equation (4.1), the divergence of equation (4.11) results in 
the pressure equation: 



The equation (4.13) discretization must use the face interpolation of H(V)/ap 



and 1/ap, which results in the following expression: 



=j;s.M) (4.14) 



4.3 OPENFOAM 

OpenFOAM (Open Source Field Operation and Manipulation) is a fiow solver of 
choice because it has a pre-existing, robust mesh motion capability that satisfies 
the GCL, and its source code is freely available through the GNU General Pub- 
lic License. OpenFOAM is a free, open source CFD software package produced 
by a commercial company, OpenCFD Ltd. It is an object-oriented library writ- 
ten in the C++ language, developing for the customized numerical solvers, and 
pre- /post-processing utilities for the solution of continuum mechanics problems. 
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including computational fluid dynamics (CFD) [2]. 



In the commercial and academic organizations, it has a large number of user 
groups across the engineering and scientiflc flelds. OpenFOAM can be applied to 
a wide range of capabilities to solve complex fluid flow, involving chemical reac- 
tions, turbulence and heat transfer. It has a set of third-party packages ParaView, 
which is used to post-process the CFD geometry and display analysis using the 
GUI [30j. OpenFOAM versions 1.6-dev has been used for this research. 



° ParaView 3.8.0 
file Edit View Sources Filters Tools Macros Heip 



1 g ^ 1 1*) pa I ^ 1 1 raii^ \Sf m w oy ^ y> w ^ :Time:[o 



ins 



3[ 



surface R : |g i-'i Si tS- fci ii [pf@| ^ 




□ Extra poi ate Pstclies 

□ Use VTKPoiyliedron 

□ Update GUI 



SMesli Parts 



_ internaiMesIn 

□ leftjniet- patch 

□ riglitjniet - patcli 

□ outiet - patcli 

□ fixedWalls-patcli 

□ frontAndBacl< - patcii 



Figure 4.4: ParaView as frontend to OpenFOAM 



1. Pre-processing: 

OpenFOAM provides a mesh generator blockMESH which the user can 
divide the geometry into many meshes. Figure (4.4) shows the geometry 
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of a square plate and we saw the use of blockMESH for mesh generation, 



shown in Figure (4.5). OpenFOAM also support converting the format of 



other CFD packages to the OpenFOAM format. 

2. Solvers: 

OpenFOAM contains solvers for incompressible flow, channel flow, com- 
bustion and stress analysis. In addition, users can create custom solvers 
without having to modify and recompile the source code with the existing 
solver. 

3. Post-processing: 

A plug-in ParaView is used to process the results of simulation cases, provide 
graphical post-processing and display analysis using the GUI. 




(a) OpenFOAM mesh (b) OpenFOAM surface with edges 

Figure 4.5: OpenFOAM block mesh 



4.4 ICOFOAM 



The flow problems to be modeled in this research are treated as incompressible and 
transient, and therefore OpenFOAMs icoFoam solver provides a good starting 
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point for the solver development. Since OpenFOAM provides a segregated algo- 



rithm to solve the coupled continuity (4.1) and momentum equations (4.1), which 



requires developing equations for each dependent variable and solved sequentially, 
an iterative method is required to solve the systems of algebraic equations. The 
icoFoam solver uses the PISO (Pressure Implicit with Splitting of Operators) 
algorithm to handle the pressure-velocity coupling. It relates to a momentum 
predictor and a correction loop, in which a pressure equation based on the volu- 
metric continuity equation is solved and the momentum is corrected based on the 
pressure change. The PISO algorithm can be described as follows: 



1. The momentum equation (4.11) is solved first by applying the estimated 
value of pressure field. Accurate source of the pressure gradient at this stage 
is unknown and the pressure field at the previous time-step is replaced. This 
stage is called the momentum predictor. The solution of the momentum 
equation gives an approximation of the new velocity field. 

2. Using the predicted velocities, the H(V) operator can be substituted and 



the pressure equation (4.14) can be evaluated. The pressure equation solu- 
tion provides the first estimate of the new pressure field. This step is known 
as pressure solution. 

3. It provides a new set of pressure field, which has always been a conservative 
fiux. As a consequence of a new pressure distribution the velocity field 
should be corrected explicitly by a velocity correction. This is the explicit 
velocity correction stage. 



The velocity corrector consists of two parts: a correction due to the change in 
the pressure gradient and the transported infiuence of corrections of neighboring 
velocities. The fact that the velocity correction is explicit means that the latter 



45 



part is neglected. The whole velocity error is assumed to come from the error in 
the pressure term. It is, however, not true and therefore is necessary to correct 
the H{V) term, formulating another pressure equation and repeating the proce- 
dure. In other words, the PISO loop consists of an implicit momentum predictor 
followed by a series of pressure solutions and explicit velocity correctors. This 
loop is repeated until the total variation in the velocity field from one time level 
to the next is less than a pre-determined tolerance. 



4.5 MYICOFOAM 

However, the fiow is treated as nonisothermal which requires solving the Navier- 
Stokes equations coupled to the energy equation. A transient solver for incom- 
pressible, laminar fiow of Newtonian fiuids myicoFoam is configured to model the 
nonisothermal reversed stagnation-point fiow in OpenFOAM. The myicoFoam 
solver is an extension of icoFoam such that it enables solving the Navier-Stokes 



equations coupled to the energy equation (4.3) 



The main framework is same as the incompressible fiow, thereby the steps is as 
follows: 



1. The momentum equation (4.11) is solved first by applying the estimated 
value of pressure field. Accurate source of the pressure gradient at this stage 
is unknown and the pressure field at the previous time-step is replaced. This 
stage is called the momentum predictor. The solution of the momentum 
equation gives an approximation of the new velocity field. 

2. The energy equation is solved in which the fiux is from solving the previous 



the momentum equation (4.11). The corresponding solution is within the 
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PISO loop, which imphes that the energy equation is solved again using the 
new flux when the new flux is evaluated. This stage is called the energy 
solution. It is speculated that the thermal coupling is as important as the 
coupling between the same pressure and velocity. 

3. Using the predicted velocities, the H(V) operator can be substituted and 



the pressure equation (4.14) can be evaluated. The pressure equation solu- 
tion provides the flrst estimate of the new pressure fleld. This step is known 
as pressure solution. 

4. It provides a new set of pressure fleld, which has always been a conservative 
flux. As a consequence of a new pressure distribution the velocity fleld 
should be corrected explicitly by a velocity correction. This is the explicit 
velocity correction stage. 



A flow diagram of myicoFoam solver is shown in Figure (4.6). The coupling of 
pressure and velocity is more important than the coupling with the density and 
temperature. Before the momentum predictor the density predictor is performed. 
Pressure solution needs the density change so the energy equation is solved flrst 
to update the density which is the main driving force for flow. 



4.6 MESH DISTRIBUTION 

While simulating the flow near the wall, due to the viscosity effect near the wall, it 
is necessary to divide much more mesh in the region near the wall. In OpenFOAM, 
the mesh distribution can be selected either uniform or non-uniform. The mesh 
deflnitions are contained in a list named blocks, consisting of a list of vertex labels, 
the number of cells in each direction and the cell expansion ratio in each direction. 
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Velocity field (divergence free) available at time n 

1 i . 

Compute the new temperature r''^^ 

I 

Compute intermediate velocities u"^ and p ' 



Solve the Poisson equation for the pressure correction p^^'""^^^^ 
u"^' is obtained from u 

I ^ 1 

Solve the velocity correction equation for u^'"^^^ ' 

u^* is obtained from u 



1 

Compute the new velocity i/''"^^and pressure/?''"^ ^fields 



Figure 4.6: Flow diagram of myicoFoam solver 
The meshes are defined as follows: 

1 convertToMet ers 1; 

2 vertices 

3 ( 

4 (0 0) 

5 (10 0) 

6 (110) 

7 (0 10) 

8 (0 1) 
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9 (10 1) 

10 (1 1 1) 

11 (0 1 1) 

12 ) ; 

13 blocks 

14 ( 

15 hex (0 1 2 3 4 5 6 7) 

16 (200 400 1) 

17 simpleGrading (1 5 1) 

18 ) ; 



The blockMesh dictionary defines a block and the mesh from the vertices, hex 
means that it is a structured hexahedral block. (0 1 2 3 4 5 6 7) is the vertices 
used to define almxlmxlm block. These sequences are very important - 
they should follow the right-hand system. (200 400 1) is the number of mesh cells 
in each direction. 



Expansion ratio = 

Axf 

Axs Axf 



Expansion Direction — ► 
Figure 4.7: Expansion ratios in a given direction 

In order to simulate the fiow near the wall, the mesh applied in the model is chosen 
to be non-uniform. In OpenFOAM, simpleGrading (1 5 1) is the expansion ratio. 
The ratio is that of the width of the final mesh Ax j along one edge of a block to 



the width of the start mesh Axs along that edge, as shown in Figure |4.7[ The 
expansion ratio allows a mesh refinement in particular direction. In our model 
the ratio of mesh widths along x— and z— axis is 1, along y— is 5. 
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Figure 4.8: Mesh distribution in the model 

4.7 CONVERGENCY 

In general, to obtain a more accurate solution, more meshes should be used in a 
numerical simulation. One should be kept in mind is that, more time is required 
to compute a solution if the domain of problem is divided into more meshes. 
In practical numerical simulations, although an accurate solution is desired, the 
number of meshes cannot be indefinitely increased because of the limitation of 
computing facilities and time constraints. To ensure the accuracy and efficiency of 
a practical numerical simulation, it is necessary to increase the number of meshes 
until no significant difference of the solutions is obtained by two consecutive sim- 
ulations. Unfortunately, this is rather difficult to estimate the optimum number 
of meshes. It takes time and patience to anticipate the mesh number. 



50 



On the other side, in numerical simulation the Courant-Friedrichs-Lewy condi- 
tion (CFL condition) is a necessary condition for convergence in solving hyper- 
bolic PDEs numerically [3lj. It is applied when explicit method is required in the 
numerical solution. When a CFD program is running, in order to achieve time 
accuracy and numerical stability, it requires that the Courant number Co in the 
flow field is always smaller than 1. The Courant number is defined for one cell as: 

Ax 

where u = 0.5 m/s is the fiow velocity in the model and Ax is the length interval. 
We therefore select based on the worst case which is the maximum At correspond- 
ing to the combined effect of a large fiow velocity and small length interval Ax. 
In our model the maximum mesh size occurs near the outlet and is equal to the 
width of the final mesh Axy along ^— axis: 

. block length (1 m)(5) 

Ax X expansion ratio = = 0.0125 m 

number oj mesh 400 

As a consequence, to achieve a Courant number less than or equal to 1 throughout 
the domain, the time step At must be less than a specified time in the time- 
marching computer simulations, otherwise the simulation will produce incorrect 
results. The time step At must be set to less than or equal to: 

CoAx (l)(0.0125m) 0.0125 

At = < ^ ^ = s 

UQ UQ UQ 

After the progress of the simulation, we hope to get accurate results in the early 

time interval, so we can later view with a post-processing package. The factor of 

0.08 is taken from the experience that can be advantageous for accuracy and thus 

in our analysis the times step is equal to 

CoAx (l)(0.0125m) 0.0125 
At = < ^ ^ = s. 

UQ UQ UQ 
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4.8 BOUNDARIES 

In this section we discuss the way in which boundaries are treated in reversed 
stagnation-point flow. In order to solve the governing equations by the numeri- 
cal method described in the previous section, boundary conditions must be pre- 
scribed. The boundaries involved in our model are not only simple geometric 
boundary conditions, but also the integral part of the solution and numerical 
simulation through boundary conditions or inter-boundary connections. 



We flrst need to consider setting up a numerical conflguration of the simula- 
tion; the boundary has to be specifled. The conditions consist of two inflow 
boundaries, an outflow boundary, and a symmetry plane on one of the two faces 
parallel to the plane of the paper and no-slip walls for the remaining boundaries. 



A schematic diagram of the problem is given in flgure 4.9, In OpenFOAM the 



boundary conditions of our problem is deflned in the hlockM eshDict dictionary: 

1 patches 

2 ( 

3 patch left.inlet ((2 6 5 1)) 

4 patch right.inlet ((0 4 7 3)) 

5 patch outlet ((3 7 6 2)) 

6 wall fixedWalls ((1 5 4 0)) 

7 empty front AndBack ((0 3 2 1) 

8 (4567)) 



9 

V 



We select the uniform velocity proflle for the inflow boundary. Discretization 
(4.14) of the momentum equation (4.10) involves the values of velocity on the 



52 



boundary. These velocity values are obtained from a discretization of the bound- 
ary conditions of the continuous problem. 

1. Inflow conditions: 

On an inflow boundary the velocities are explicitly given; we impose this for 
the velocities normal to the boundary by directly fixing the values on the 
boundary line. 

2. Outflow conditions: 

In the outflow boundary condition the normal derivatives of both velocity 
components are set to zero at the boundary, which means that the total 
velocity does not change in the direction normal to the boundary, i.e., 

du dv ^ 
dx dy 

3. No-slip condition: 

The continuous velocities should vanish at the wall boundary to satisfy the 
no-slip condition. For the values laying directly on the wall boundary we 
thus set both velocity component to zero. 

u{x^ 0) = v{x^ 0) = 

4. Symmetry plane: 

Our problem is a two-dimensional problem which is symmetric about z— axis. 
This means that boundary condition refers to a planar boundary surface. 
Values lying directly on the boundary are not required to calculate. 

The expression of the velocity values on the boundary is shown as following: 
Listing 4.1: Velocity Boundary conditions 

1 boundaryField 
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2 { 



10 



11 



left_iiilet {"type 
value 

right_iiilet {type 
value 

outlet {type 

fixedWalls {type 
value 

f ront AndBack {type 



f ixedValue ; 
uniform (-1 0);} 
f ixedValue ; 
uniform (10 0);} 
zeroGradient ; } 
f ixedValue ; 
uniform (0 0);} 
empty;} 




^^7/ / // // // // // // // // // // // / / / 



Figure 4.9: Boundaries of reversed stagnation-point flow in a control volume 



Each patch deflnes a type, a name, and a list of boundary faces. The patch is 
defined by three sides of the block based on the vertex numbers. The order of the 
vertex numbers is such that they are marched clockwise when looking inside from 
the control volume. For example, type of the boundary fixedWalls is defined as 
wall. For two-dimensional fiow the fiow field on the boundary front AndBack is 
not required to evaluate, thereby is defined as empty. 



For the temperature profile in nonisothermal fiow, there are essentially two dif- 
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ferent boundary conditions; to impose these, we divide the boundary into two 



1. Dirichlet boundary conditions: 

Using this boundary condition, the constant wall temperature is pre- 
scribed at the wall. The temperature of the fluid from a wall may be de- 
scribed in the form 



2. Neumann boundary conditions: 

This boundary condition describes how much heat is passed on to the wall 
by the fluid. This is determined by both the material properties of the 
wall and the temperature difference across the wall. For a constant fluid's 
thermal conductivity k and heat flux q^j across the wall, it may be described 
in the form 



In our model, the external flow temperature and wall temperature are constants 
and Dirichlet boundary conditions are required. The expression of the tempera- 
ture values on the boundary is shown as following: 



parts: 



T{x,0)^T, 



W 




Listing 4.2: Temperature Boundary conditions 



1 boundaryField 



2 boundaryField 



3 { 



4 



lef t_inlet 



{type 



f ixedValue ; 



5 



value 



uniform 373;} 



6 



r ight_inlet 



{type 



f ixedValue ; 



7 



value 



uniform 373;} 
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8 outlet {"type zeroGradient ; } 

9 fixedWalls {"type fixedValue; 

10 value uniform 273;} 

11 front AndBack {type empty;} 

12 } 

V J 



On the other side, the selection of fluid is diflicult because crude oil is predomi- 
nantly a mixture of hydrocarbons. Under surface pressure and temperature con- 
ditions, the lighter hydrocarbons methane, ethane, propane and butane occur as 
gases, while the heavier ones from pentane and up are in the form of liquids or 
solids. 



It is, however, in the underground oil reservoir the proportion which is gas or 
liquid varies depending on the subsurface conditions. This represents that the 
flow system is not single-phase, but is multiphase. As a result, in order to sim- 
plify the diflSculties in simulation, it is possible to select another fluid to replace 
the crude oil. After simulation, the numerical result can be analyzed in the crude 
oil situation, by comparing Reynolds number. We create a fluid that has sim- 



ilar physical properties of crude oil. In Fig. (4.10), the viscosity of crude oil is 
approximately 10 cp , which equals to 0.01 Pa • s. 



4.9 DIMENSIONESS NUMBER 

Besides, it is required to know fluid physical properties. For our transient solver 
myicoFoam the physical properties nu and DT are stored in the transport Properties 
flle which is a dictionary for the dimensioned scalar. The flrst items loaded is the 
kinematic viscosity from the transport Properties dictionary flle and is equal to 
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Figure 4.10: Properties of crude oil in enhanced recovery methods [3] 



u. Another transport property related to the thermal diffusion denoted as DT 
equals to DT = u/Pr 



In this thesis, effect of reversed stagnation-point on the nonisothermal flow fleld 
behavior is studied. The two parameters uq and DT^ which are the inflow velocity 
at both the left and right boundary and molecular thermal diflFusivity respectively, 
are investigated in the numerical simulation. 



One of the contributions in this thesis is to acquire the relationship of Reynolds 
number and flow velocity, where the reversed stagnation-point flow exists under 
the condition of ensuring the flow is laminar. Reynolds number is a dimensionless 
flow property. It gives a measure of the ratio of inertial forces to viscous forces 
and, consequently, it quantifles the relative importance of these two types of forces 
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for given flow conditions. Reynolds number is defined as: 

V 

where V is the mean velocity and L is the characteristic length, equals to half of 
the length of wall. Reynolds number can also describe the property of the flow, 
whether it is laminar, transition or turbulent flow. For a smooth flat plate with 
a uniform free stream, the transition process begins at a critical Reynolds num- 
ber, Rccritical ~ 1 X 10^, and continues until to the turbulent at the transition 
Reynolds number, Rctransition ~ 3 x 10^. The flow is said to be laminar flow 
when Re < RCcritical ~ 1 x 10^. Several cases of simulations with varied Reynolds 
number were performed, in which the value of i/q is chosen from 0.5 to 20 with 
flxed value of v in part of the simulations. 



Reynolds number can be obtained as applying the nondimensional form of the 
incompressible Navier- Stokes equations: 

' dV 



P 



dt 



+ V)V 



= -Vp + /iV^V + / 



(4.15) 



When the equations undergo the dimensionless analysis, that is when it is mul- 
tiplied by a factor with inverse units of the origin equation, we acquire a form 
which does not depend directly on the physical sizes. One possible way to get a 
nondimensional equation is to multiply the whole equation by the factor L/pV"^ 
and to set 

The Navier-Stokes equation can be rewritten without dimensions: 



+ {V' ■ V)V' = -VV + V'2 W + f 



(4.16) 
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Finally, dropping the primes, we have 

^ + (V • V)y = -Vp + ^V^F + / (4.17) 

Therefore, all flows with the same Reynolds number are comparable mathemat- 
ically. It is noted that, in the above equation, as Re ^ oc the viscous terms 
vanish. High Reynolds number flows are approximately inviscid in the external 
flow. Meanwhile, the velocity components and nondimensional variable of our 
problem can be rewritten in the form: 



u 



-xfM (4.18a) 



V = Re" 2/(7/) (4.18b) 
T] = Re^y (4.18c) 
r^t (4.18d) 



and the governing similarity equation remains unchanged: 

fr]T — ifr]) + ffrjT] — frjrjr] — (4.19) 

That is why the results of direct numerical simulations and that of similarity anal- 
ysis are comparable with the same Reynolds number. On the other side, studying 
the effect between different Prandtl number and temperature distribution in non- 
isothermal flow is the second goal in this thesis. In Chapter 2 the deflnition of 
Prandtl number is introduced as: 



k DT 

Prandtl number is the ratio of momentum diffusivity to thermal diffusivity. Sev- 
eral cases of simulations with different Prandtl number were performed, ranging 
from 0.003 to 0.1 with flxed value of u in part of the simulations. The following 
table illustrates the parameters used in this study. 
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Table 4.1: Reversed stagnation-point flow with parameters 

Case y DT uq (m/s) At (s) Re Pr 

1 0.01 0.01 1 0.001 50 1 

2 0.01 0.01 2 0.0005 100 1 

3 0.01 0.01 5 0.0001 250 1 

4 0.01 0.01 10 0.00005 500 1 

5 0.01 0.01 20 0.000025 1000 1 

6 0.01 0.01 50 0.00001 2500 1 

7 0.01 0.01 100 0.000005 5000 1 

8 0.01 0.01 200 0.0000025 10000 1 

9 0.01 0.033333333 20 0.000025 1000 0.3 

10 0.01 0.014285714 20 0.000025 1000 0.7 

11 0.01 0.003333333 20 0.000025 1000 3 

12 0.01 0.001428571 20 0.000025 1000 7 

13 0.01 0.001 20 0.000025 1000 10 



CHAPTER 5: RESULT AND DISCUSSION 



In this chapter, the numerical results of reversed stagnation-point flow in Open- 
FOAM will be discussed. The results of direct numerical simulations are compared 
to the analytical solutions of the reversed stagnation-point flow to ensure valida- 
tion of modeling in the simulations and to check the reliability of the numerical 
results. 



5.1 FLOW VISUALIZATION 

We can plot the position of each particle in our simulation inside of the control 
volume to see the effects of the streamlines for various Reynolds number. The 



following pages (Figures 5.1 to 5.7) show the stream lines, both evolving in time 
as well as at steady state, at various uq. At t = 0, the inflow velocity is instanta- 
neously set from zero to uq^ thereby slowly setting in motion the fluid initially at 
rest. 



When Reynolds number is relatively small, say Re < 50, convective forces can 
be neglected as compared to viscous forces and the laminar boundary layer sepa- 
rated from the wall at the reversed stagnation point. With an increase in Reynolds 
number, both convective forces and viscous forces are in the same order. The lam- 
inar boundary layer starts separating from the wall before the reversed stagnation 
point. At the same time, there emerges a symmetrical pair of stable vortices which 
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0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



(c) t = 0.2 (d) t = 0.3 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



(e) t = 0.4 (f) Steady 

Figure 5.1: Stream Line, time evolution at Re = 50 
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(a) t = 0.01 (b) t = 0.1 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



(e) t = 0.25 (f) t = 0.28 

Figure 5.2: Stream Line, time evolution at Re = 100 
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Figure 5.2: Stream Line, time evolution at Re = 100 
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Figure 5.3: Stream Line, time evolution at Re = 250 
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Figure 5.3: Stream Line, time evolution at Re = 250 
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(e) t = 0.05 (f) t = 0.055 

Figure 5.4: Stream Line, time evolution at Re = 500 
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Figure 5.4: Stream Line, time evolution at Re = 500 
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Figure 5.5: Stream Line, time evolution at Re = 1000 
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(k) t = 0.045 (1) Steady 

Figure 5.5: Stream Line, time evolution at Re = 1000 
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Figure 5.6: Stream Line, time evolution at Re = 2500 
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Figure 5.6: Stream Line, time evolution at Re = 2500 
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Figure 5.7: Stream Line, time evolution at Re = 5000 
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Figure 5.7: Stream Line, time evolution at Re = 5000 
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create a back flow, and hence, a circulation region forms close to the reversed stag- 
nation point. With a further increase in Reynolds number, the laminar boundary 
layer becomes thicker and the vortices extend. The corresponding steady state 
velocity proflle indicates that near the wall most of the fluid has a reversal di- 
rection, allowing two steady symmetric eddies to form in the resulting gap. In 
all cases, however, the flow reaches steady state, hence streamlines coincide with 
streaklines. 



The origins of boundary-layer separation are associated with the frictional forces 
within the boundary layer and a positive or adverse pressure gradient occurs in 
the direction of flow. Near the wall region, some fluid energy is dissipated in 
overcoming friction in the boundary layer. When vortices are formed on the de- 
celerated boundary layer, the flow tries to decelerate in a short manner. In the 
entire boundary layer, once the outer flow is accelerated by a pressures drop, the 
fluid elements will also move in the direction of motion, and hence, the flow will 
keep in its original direction along the surface. On the other hand, if the pressure 
of particles declines in the direction opposite to the flow, the outer flow is there- 
fore decelerated. The remaining energy is not suflicient to overcome the increased 
pressure. Then slower fluid particles of the boundary layer are even more slowed 
down. Eventually, if the deceleration is large enough such that the flow particles 
stop in motion and start moving in the opposite direction, the flow separates from 
the wall and a backflow region emerges. 



Fig. (5.8) demonstrates the vortex formation in the pressure distribution p. When 
the streamline portrait of the boundary-layer flow is close to the separation po- 
sition A, since the backflow is close to the wall, the separation rolls up into one 
or more vortices. Soon after, a great thickening of the boundary layer exists near 
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Figure 5.8: Separation process (maximum velocity M, separation point A) [4J 



this region. 



At the separation point the waU streamhne departs the waU at a certain an- 
gle. The position of the point of separation is that point on the wall where the 
velocity gradient perpendicular to the wall vanishes. In another words, the point 
where the walls shear stress becomes zero. 



5.2 VELOCITY PROFILE 



Next we discuss the velocity field frj. Figures (5.10) to (5.13) show the similarity 
velocity distribution along the ry-direction at locations ofx = 0, 0.1, 0.2 and 
0.3 respectively, evolving in time at various Reynolds number. Since the case of 
opposing fiow can be mapped to a case of fiow by x ^ — x, we will not present 
here values of the negative value of x for the case of opposing fiow. 



It can be observed from these figures that at the beginning of fiuid motion, the 
minimum value of frj almost keeps to be zero when the value of Reynolds number 




(e) X = 0.1 
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ranging from 50 to 2500. No back flow is observed near the wall region. One of 
the reasons of this phenomenon is that, at the beginning of the fluid motion, in 
just a very short period of time, the viscous forces have propagated mostly into 
the fluid. As the fluid is at rest initially, the fluid flow has to overcome a large 
inertia, resulting in a fluid flow motion. It seems that the convective forces can 
be neglected as compared to viscous forces. 

We examine two phenomena here: the dependence of flow velocity on x and the 
dependence of the external flow. One of the assumption in the analytical solu- 
tion is that the velocity fleld fj^ is a function of vj only in the region near the 
reversed stagnation point, provided that the velocity fleld f^q is independent of x. 
When the fluid flows near the origin or the location of x-coordinate is relatively 
small, say x < 0.1, the distribution of frq is independent of x. On the other side, 
the numerical solutions show variation of velocity along the x-direction. Large 
discrepancy occurs as a larger value of x is applied in the numerical simulation, 
which violate the assumption of no variation of velocity along the x-direction in 
the region far away from the reversed stagnation point. 



-t=0.01 
-t=0.03 
-t=0.05 
-t=0.07 
-t=0.09 



1 



5 



10 



(a) X = 





-t=0.01 
-t=0.03 
-t=0.05 
-t=0.07 
-t=0.09 



5 



10 



(b) X = 0.1 

Figure 5.10: Similarity velocity field, time evolution at Re = 100 
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Figure 5.10: Similarity velocity field, time evolution at Re = 100 
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Figure 5.11: Similarity velocity field, time evolution at Re = 250 
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Figure 5.11: Similarity velocity field, time evolution at Re = 250 
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Figure 5.12: Similarity velocity field, time evolution att Re = 500 
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Figure 5.12: Similarity velocity field, time evolution at Re = 500 
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Figure 5.13: Similarity velocity field, time evolution at Re = 1000 
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Figure 5.13: Similarity velocity field, time evolution at Re = 1000 
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Figures (5.14) to (5.17) show comparisons between the numerical simulations and 
the similarity solutions of reversed stagnation-point flow for different values of r. 
Lines without markers denote results obtained from numerical simulation (NS), 
dotted lines are obtained from the flnite-difference formulations (SS). In the region 
near the reversed stagnation point the solution agrees remarkably well for smaller 
values of r with the known similarity solution, thus conflrming the predictions of 
the viscous Proudman- Johnson solution. On the other hand, discrepancy occurs 
as a larger value of r is applied in the numerical simulation. However, far away 
from the wall region, a large diflFerence is observed from the results obtained by 
these two solutions. The component of velocity normal to a wall is not outward 
the wall in the region near the reversed stagnation point. The vorticity created at 
the wall will be convected outward the wall, which spreads the vorticity towards 
its source at the boundary. 



The other consideration is the behavior of external flow. Proudman and Johnson 
considered a constant potential flow outside the boundary layer. From Figures 
(5.10) to (5.17), one may observe that as ^ oc, the similarity velocity frj cannot 
ultimately approach to 1. It is clearly observed that frj gradually drops when the 



time step increases. As mentioned in previous section, near the wall region or in 
the boundary layer the phenomenon of reversed flow with boundary-layer separa- 
tion occurred. There is no justiflcation whatever outside the boundary layer for 
supposing that at large distances from the wall {rj oc) the velocity '^(x,7y,r) 
should pass over smoothly into that for inviscid Vq. Also we proofed that frj can- 
not ultimately approach to 1 from above or below, nor in an oscillatory manner 



so that no solution to equation (3.8) exists in two-dimensional steady case. As a 
consequence the assumption that the potential flow Vq is restricted not to be a 
constant as well as a time dependent function is reasonable in reversed stagnation- 
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Figure 5.14: Comparison between the numerical velocity profiles and similarity 
velocity at Re = 100 




Figure 5.15: Comparison between the numerical velocity profiles and similarity 
velocity at Re = 250 
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Figure 5.17: Comparison between the numerical velocity profiles and similarity 
velocity at Re = 1000 
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point flow. 

Moreover, two opposed vortices emerge in the regions are usually in the vicinity 
of the boundary of the fluid adjacent to wall where viscous forces are dominant. 
The most important implication of the solution contemplated is the growth of 
vortices near the wall in a main stream. According to the present similarity 
solution, separation of the main flow cannot start at any flnite time r in the limit 
as ^ 0. Moreover, the inviscid Proudman- Johnson solution implies a steady 
flow ^ — 1 when r ^ oc and the flow problem becomes the classic stagnation- 
point problem (Hiemenz [5j) by changing the sign in /. The solution shown in 



Figure (3.2) indicates that the region of reversed flow expands and has inflnite 
dimensions as r ^ oc, which violates the results of numerical simulations that 
two flnite-dimensional vortices appears near the wall in steady state. Proudman- 
Johnson solution is only approximate but cannot guarantee that it is free from 
an inflnite multiplicative error for large times. 

5.3 TEMPERATURE PROFILE 

Now we return to the numerical results of the nonisothermal stagnation-point 



flow problem. The following pages (Figures 5.18 to 5.23) show the heatlines, both 
evolving in time, at various Prandtl number. The thermal color was illustrated 
in the rainbow scale. Colors closer to red are hot areas and colors closer to blue 
are cold areas. At t = 0, the inflow velocity is instantaneously set from zero to 
uq^ thereby slowly setting in motion the isothermal fluid initially at rest. The 
heatlines in these flgures show heat flowing mainly from the cooled wall to the 
heated external flow by conduction in the beginning. The heated external flow 
passes though the wall and rises, and as it does, it cools down by conduction 
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and convection of heat. After closing to the reversed stagnation point, under the 
motion of backflow, it sinks to the wall where it is prohibited from sinking fur- 
ther. This hot fluid has thermally contracted to become dense near the reversed 
stagnation point along the edges of the wall. It trapped in the region near the 
cooled wall starts to cool down. 

It is worth talking into consideration that for liquid metals the Prandtl number 
is very small (Pr << 1, generally in the range from 0.01 to 0.001. They have a 
high thermal conductivity and low viscosity. The value of Pr = 1 corresponds to 
diatomic gases, including air. For many fluids, including water, Prandtl number 
lies in the range from 1 tolO. Large values of Pr (>> 1) correspond to high- 
viscosity oils and Pr = 7 corresponds to liquid water at room temperature. 
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Next we discuss the temperature distribution B. Figures (5.24) show the simi- 
larity temperature distribution along the ry-direction, evolving in time at various 
Prandtl number. As anticipated, since the temperature distribution is indepen- 



dent to X, we do not discuss the region where x < 0.1. From Figures (5.24), 
the temperature distribution is monotonically increasing. B drops from a remote 
value to its value inside the thermal boundary layer adjacent to the wall. Near 
the backflow region, surprisingly, no discernible temperature signature appears 
between the dividing streamlines. The dimensionless temperature B is linear pro- 
portional to the dimensionless distance rj in the region close to the wall, which is 
consistent to our similarity temperature solution. 



It is noticed that the dimensionless wall temperature gradient B^(0) raises with 
increase of Prandtl number, but the thermal boundary layer thickness decrease 
with increase of Prandtl number. Larger Prandtl numbers results in the thin- 
ner boundary layers and larger temperature gradients near the wall. When Pr is 
small, the heat diffuses very quickly compared to the velocity field and hence for 
liquid metals the thickness of the thermal boundary layer is much thicker than 
that of the velocity boundary layer. 



We are interested in comparison between the numerical simulation and the similar- 



ity result. Figures (5.25) to (5.25) illustrate comparisons between the numerical 
simulations and the similarity solutions of nonisothermal reversed stagnation- 
point fiow, when the dimensionless Reynolds number Re = 1000. Lines without 
markers denote results obtained from numerical simulation (NS), dotted lines 
are obtained from the finite-diflFerence formulations (SS). It is shown that when 
Prandtl number is less than 1, our simulated results fall within the values obtained 
from the finite-diflFerence formulations. 
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Figure 5.24: Similarity temperature field, time evolution 
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Figure 5.24: Similarity temperature field, time evolution 
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Figure 5.24: Similarity temperature field, time evolution 
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Figure 5.25: Comparison between the numerical temperature profiles and simi- 
larity temperature 
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Figure 5.25: Comparison between the numerical temperature profiles and simi- 
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CHAPTER 6: PARTICULAR SOLUTION 



We complete the discussion of the Proudman- Johnson equation. Our objective is 
to obtain a similarity solution of the governing equation. Comparing to the results 
of numerical simulation, it is found that the potential flow Vq may be expressed 
as a time dependent function. In this chapter, rather than considering inviscid 
flow as external flow, it could instead be thought of a monotonic potential flow 
in balance of both viscous and convection terms in the total flow fleld. Let us 
discuss the similarity solution in a diflFerent manner. 

6.1 ANALYTICAL ANALYSIS 

In our two-dimensional model, the fluid remains at rest when time t < and is set 
in motion at t > such that at large distances far above the planar boundary the 
potential flow is a constant Vq for all value of t. Both Proudman and Johnson [T], 
and Robins and Howarth [H] have set Vq = 1 and the corresponding boundary 
condition /7^(oc,r) = 1. When the flow is in steady state such that frjr = 0, it 
was proven that the similarity velocity frjiv) cannot ultimately approach to 1. 
The differential equation has no solution. Smith [15j generalized the solution of 
Proudman and Johnson with both viscous and convection terms in balance by 
considering the monotonic potential flow Vq ^ot to be a constant when the time 
is relatively large. 
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Ill 



When the flow decays so rapidly that viscous force cannot be ignored away from 
the boundary, the viscous terms term frjrjrj must be included in the entire flow 
field. If the potential flow Vq is restricted not to be a constant, then the boundary 
condition /7^(oo,r) may be expressed in a time dependent function. Numerical 
solution of reversed stagnation-point flow for this particular case has been studied 
in [24]. 



Now we go through the analysis of this particular case. As with the governing 
equation of reversed stagnation-point flow, we can write the stream function as 

ip = —^/Auxf{ri^T) (6.1a) 

77 = Y ~y (6.1b) 
T = At (6.1c) 

where ^4 is a constant proportional to Vo(r)/L, Vo(r) is the external flow velocity 
removing from the plane and L is the characteristic length. These result in the 



governing equation (2.15) 



[fr]T - {fr]f + ffrjT] - frjrjr]] = function of r only, 
or the function of r may be expressed as 

fvr - Unf + ffm - fvm = -C{r), (6.2) 

Under the boundary conditions /7^(oo,r) = 1, the value of C{r) should be a 
constant and equal to 1. If the boundary condition /7^(oc,r) is restricted not 
to be a constant, following the assumption of Shapiro |[16j, a particular time- 
dependence function C{r) may be expressed in the form 

C{r) = 4 (6.3) 
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where c is an arbitrary constant. The partial differential equation can be simpli- 
fied by a similarity transformation when a new similarity variable is introduced. 
This converts the original partial differential equation into an ordinary differential 
equation. 



When T is small the solution may be obtained by the method developed by Bla- 
sius [32j and the solution satisfying the early stages of the diffusion are of the 
form 

/ = X function of ( ^ ) (6.4) 



For small values of r, therefore, the variable rij^Jr is more appropriate than 77 



itself. When we consider equation (6.2), if a time dependent function is taken 



into account, the diffusion variable transformation is introduced 

^ = ^ (6.5a) 
/(ry,r) = ^F(0 (6.5b) 

Here ^ is the time combined nondimensional variable and F is the nondimensional 
velocity function; F is then the sole function of ^ and insertion of the similarity 
transformation yields an ordinary differential equation 

- \f^^ -F^ - F'^ + FF" - F'" = -c (6.6) 

Li 

where the prime denotes the derivative with respect to the variable ^. 



Equation (6.6) is a third-order nonlinear ordinary differential equation. A crucial 
step in obtaining an analytical solution involves rearranging the equation as an 
autonomous differential equation. In mathematics, an autonomous differential 
equation is a system of ordinary differential equations which does not explicitly 
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depend on the independent variable. 



In order to omit the variable ^ in the differential equation, it is generally accepted 
as a change of variable 



and the equation becomes to an autonomous differential equation 

3 



(6.7) 



QQ" - 2Q' - g'2 _ Q'" = _c + 



(6i 



In our analysis, P = Q' is the dependent variable and Q is the independent 



variable. Equation (6.8) is reversed ranged as 



QP' - 2P - p2 _ p" = -c + 



(6.9) 



and the chain rule reduces equation (6.9) to a second-order ordinary differential 
equation 



Equation (6.10) is analytically solvable that the solution might be expressed as a 



low order polynomial. It is suggested that 



p = a + bQ + dQ' 



(6.11) 



and substituting into equation (6.9) and comparing the coefficients in the powers 



of Q results in a system of linear algebraic equations 

2a^d + ab^ + 2a-a^ = -c + 
8abd + b^ + 2b + ab = 
8ad'^ + ilb^ + 2)d = 
I2bd^ -bd = 
6d^ -d^ = 



(6.12a) 
(6.12b) 
(6.12c) 
(6.12d) 
(6.12e) 
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Solving the related algebraic equation, we have 



a = —- 



b = 0, c 



d 



6 



(6.13) 



Substituting the constant into equation (6.11) yields a first-order differential equa- 
tion 

(6.14) 



2 6 



Equation (6.14) is Riccati equation, which is any ordinary differential equation 



that is quadratic in the unknown function. The standard form of classical Riccati 
equation is 

Q' = RQ^ + SQ + T (6.15) 

The solution of Riccati equation can be obtained by a change of dependent vari- 
able, where the dependent variable y is converted to q by [33] 



Q 



q R 



(6.16) 



1 ^ 

By identifying i? = = and T = — 2, the change of variables in equation 



(6.14) becomes 



6(7' 



(6.17) 



SO the equation (6.14) becomes a second-order linear differential equation 



(6.18) 



of which the general solution is 



q = A cosh - + B sinh - 



(6.19) 



where A and B are arbitrary constants. Applying this solution in equation (6.7) 



leads to the general solution of equation (6.6) 



(; 3 A sinh ^ + 3B cosh ^ 
2 A cosh 1 + 5 sinh | 



(6.20) 
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Application of the impermeabihty condition F{0) = leads to the determination 
of the constant 5 = 0, so the exact solution becomes 

F(0 = |-3tanh| (6.21) 

Collecting results, the velocity function becomes 

/(^,r) = ^(|-3tanh|) (6.22) 



/ A 

where <^ = \ — y is the non-dimensional distance from the plate. In view of 
V lyr 



(6.22), the flow far away from the boundary becomes 



lim firj, = ^ - 3) = f - 4 (6.23) 



where tends exponentially to a positive constant as <^ ^ oc. The flow fleld is 
not able to remain unchanged at suflicient distances far away from the wall at any 
finite time, the potential flow cannot be assumed as the outer boundary condition 
for all values of r. A continuous change as Vo(r) decreases in magnitude for large 
r should be expected outside the boundary. 



Our objective is to obtain a particular solution of the unsteady reversed stagnation- 
point flow. The solution is obtained in the similarity transformation for unsteady 



viscous flows. The flrst term of (6.23) shows that the external flow is directed 
toward the ^— axis and away from the wall. The appearance of a negative value 



in the second term in (6.23) describes a uniform velocity directed toward the wall. 
The function tanh<^ has a Taylor series expansion with only odd exponents for ^, 
that is 

2c^ 17c^ 

tanlK = .-3+3^-^ + ... (6.24) 
Thus, the flow near the boundary becomes 



lini /(?7,r) 



< 



r 



(6.25) 
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Surprisingly, the component of velocity normal to a wall is not outward the wall 
in the region near the reversed stagnation point. The vorticity created at the 
wall will be convected outward the wall, which spreads the vorticity towards its 
source at the boundary. An explanation is that an adverse pressure gradient in 
the region close to the wall leads to a boundary-layer separation and associated 
flow reversal, and therefore the flow divides into a wall region of reversed flow and 
an outer region of forward flow. 



At this part it is particular to emphasize a point which seems to been ignored 
in the analysis. Near the wall region or in the boundary layer the phenomenon 
of reversed flow with boundary-layer separation occurred. No trouble arose from 
the idealization of Proudman and Johnson that the viscous forces are of the 
same order as the inertial forces near the stagnation point. Since no information 
concerning the nature of the flow for flnite times has yet been included, there is 
no justiflcation, theoretical or experimental, for supposing that at large distances 
from the wall {t] oc) the velocity v^x^rj^r) should pass over smoothly into 
that for inviscid Vq. Once the reversed flow has occurred, the external boundary 
condition must be affected and that the whole problem becomes conceptually 
unsound. 



6.2 NUMERICAL SOLUTION 
6.2.1 VELOCITY DISTRIBUTION 



The particular solution (6.22) is noteworthy in that it is completely analytical. 
Now this solution satisfles the Navier-Stokes equations; however the equation has 
no solution that satisfles the necessary no-slip condition at the wall in the pres- 
ence of non-zero term F\0) = —1. 



117 



In order to satisfy this too, the effect of no-shp condition F\0) = must be 
taken into account. To do this we apply the numerical analysis for the velocity 
distribution. The similarity equation and the relevant boundary conditions are 

-1<;F" -F' - F'2 + FF" - F'" = -c 

F(0) = F'(0) = (6.26) 
F'(oo) = 1 

where c = 3/4 to satisfy the unsteady viscous flows in the outer region. 



Equation (6.26) is a third-order nonlinear ordinary diflFerential equation. It is 



convenient to describe the problem in terms of a system of flrst-order equations 
when solving an ODE system numerically [34]. In numerical analysis, the Runge- 
Kutta methods are an important family of implicit and explicit iterative methods 
for the approximation of solutions of ordinary differential equations. This method 
applies a trial step at the midpoint of an interval to cancel out lower-order error 
terms, besides; Runge-Kutta formulas are the methods of solving initial value 



problems for ordinary differential equations. Since (6.26) is a boundary- value 



problem, apparently we have to alter the boundary value conditions into the ini- 
tial value conditions. 



For example solving an n^^-order problem numerically is common practice to 
reduce the equation to a system of n flrst-order equations. Then, by deflning 
yi = F, y2 = F' ^ Us = F" ■> the ODE reduces to the form 



y2 



(6.27) 
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The first task is to reduce the equation above to a system of first-order equa- 
tions and define in MATLAB a function to return these. The relevant MATLAB 



expression for equation (6.27) would be: 



Listing 6.1: System of first-order equations 



1 function dy = stagnat ion (t , y ) 

2 c = 3/4; 

3 dy = zeros (3,1) ; 
dy(l) = y(2); 
dy(2) = y(3) ; 

6 dy(3) = c-l/2*t*y(3)-y(2)-y(2)*y(2)+y(l)*y(3) ; 



end 



The next step is to convert the boundary value into initial value, because orfe45, 
an ODE solver in MATLAB, can only solve the initial- value problem. From 



equation (6.26), we gauss the value of F^\0) such that F\oo) = ^. The commands 
written in MATLAB would be 

Listing 6.2: ODE solver 



1 function main 

2 [T,Y] = ode45(@stagnation , [0 10], [0 

3 end 



•1]); 



The complete solutions of two-dimensional stagnation-point fiow with diflFerent 
values of F^\0) are shown from Figures (6.1) to (6.1). In these figures the sim- 



ilarity stream function F, the velocity profile and the shear stress F^^ are 
represented. This solution is a similarity solution of the reversed stagnation- 
point fiow over a fiat plate, describing an unsteady viscous fiow in both outer and 




(b) F"(0) = 

Figure 6.1: Numerical solutions of viscous reversed stagnation-point flow 




(d) F"(0) = -1 

Figure 6.1: Numerical solutions of viscous reversed stagnation-point flow 




(f) F"(0) = -1.7 

Figure 6.1: Numerical solutions of viscous reversed stagnation-point flow 
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inner regions. 



The result looks interesting from both theoretical and engineering points of view. 
A single dividing streamline plane separates streamlines approaching the plate 
from external flow streamlines. The boundary-layer thickness increases as the 
square root of r. The boundary layer thickness is the distance from the body at 
which the velocity is 99% of the velocity obtained from an inviscid solution. When 
F^\0) > 0, the values of F and are always greater than zero. No separation 
occurs near the wall region. 



The similarity velocity flelds are shown in Figures (6.2) at diflFerent values of 
F^\0). It is reasonable to state that, in general, separation will occur near the 
wall as ^ and the region of reversed flow will move outward away from the 



wall as F^\0) < —1. Moreover, it is noted that given from equation (6.23) the 
external flow velocity 

3 

fiv^r) = 

2t ^Jt 



will tend to zero for large times r ^ oc. We have, from equation (6.26) with Vq 
= 0, the equation 



F(0) = F'(0) = 
F'(oc) = 



(6.28) 



where c = 0. Figure (6.3) shows the numerical solutions at various values of 
F^^(O), indicating that the nonlinear convective terms play a secondary role in 
fluid motion as r ^ oc, the viscous forces may play a signiflcant role to decelerate 
the velocities to zero. The boundary of this region comes to rest and flnally the 
region of reversed flow does not continue to grow but has flnite dimensions. Larger 
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(b) F"(0) = 

Figure 6.2: Similarity velocity field as a function of <j 
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(d) F''(0) = -1 
Similarity velocity field as a function of ^ 
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Figure 6.2: 



(f) F"(0) = -1.7 
Similarity velocity field as a function of <j 
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value of F^\0) corresponds to larger dimension of the reversed region. 



6.2.2 TEMPERATURE DISTRIBUTION 



Under the assumption that the viscous dissipation is negligible compared to con- 
duction at the wall, 9 = Q{(;) is the function of (; only. The energy equation may 
be written as 

e'' + Pr - = (6.29) 

subject to the boundary conditions 



9(0) = 0, e(oc) = 1 



(6.30) 



where Pr = — is the Prandtl number. Equation (6.29) is a second-order linear 
k 



ordinary differential equation, and has an exact solution through a transformation. 
Let 



Z = 



de 



(6.31) 



Substituting equation (6.31) into equation (6.29) and simplifying gives 



dZ _ 
d(; 

A further integration provides 

Z = Zq exp 



■Pr[-i-F]Z 



or 



Q = Zq exp 



-Pr (^^ ~ 1 



+ eo 



Compare to the boundary conditions, we get 



1 



oo 



eo = o 

d(^ exp —Pr j ^-5 — F \ ds 




Figure 6.3: Numerical solutions of viscous reversed stagnation-point flow for c = 
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An exact solution of equation ( |6.29 ) is given as 

fn exp —Pr fn ( is — F I ds 

= ^ T ^^-22) 

Jq^ rf^ exp —Pr Jq ~ Fj ds 

A closed-form solution of the thermal energy equation for forced convection sys- 
tem is obtained. The solution, however, is not anticipated to integrate because 



equation (6.21) does not satisfy impermeability condition of the wall = and 



we cannot have an analytical solution of F. It is convenient to solve the decou- 
pled momentum and energy equations numerically. Defining ?/4 = B, 7/5 = 6^ 



and combining the variables in the momentum equation (6.27), the uncoupled 
momentum and energy equations reduce to the form 



d<. 



2/2 
?/3 



c - 5^?/3 -y2-yl + yivs 



(6.33) 



The relevant MATLAB expression for (6.33) would be: 



Listing 6.3: System of first-order equations 



1 function dy = stagnation (t , y , Pr) 

2 c = 3/4; 

3 dy = zeros (5,1) ; 
dy(l) = y(2); 
dy(2) = y(3); 

6 dy(3) = c-l/2*t*y(3)-y(2)-y(2)*y(2)+y(l)*y(3) ; 

dy(4) = y(5) ; 
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8 dy(5) = Pr*y(5)*(y(l)-l/2*t) ; 

9 end 



As was previously indicated, the boundary value problem is changed into initial 
value problem by taking a gauss of 0^^(O) such that B^(oc) = 1. The corre- 
sponding commands written in MATLAB would be 

Listing 6.4: ODE solver when Pr is equal to 1 

1 function main 

2 [T,Y] = ode45(@stagnation , [0 10], [0 -1.03 0.8]); 

3 end 



The numerical solution for temperature distributions is shown in Figure (6.4). 
It is noticed that the dimensionless wall temperature gradient B^(0) raises with 
increase of Prandtl number, but the thermal boundary layer thickness decrease 
with increase of Prandtl number. The thermal boundary layer thickness is the 
distance from the body at which the temperature is 99% of the temperature 
obtained from an inviscid solution. The decrease of thickness can be explained 
by the definition of Prandtl number that is inversely proportional to the thermal 
diffusivity k/pcpu. If the Prandtl number is greater than 1, the thermal boundary 
layer is thinner than the velocity boundary layer. If the Prandtl number is less 
than 1, which is the case for air at standard conditions, the thermal boundary 
layer is thicker than the velocity boundary layer. 



130 




Figure 6.4: Reversed stagnation-point temperature distributions B for various 
value of Prandtl number 



CHAPTER 7: CONCLUSION AND 
RECOMMENDATION 



In this study, nonisothermal stagnation-point flow is studied by applying an un- 
steady numerical model in Computational Fluid Dynamics. Beyond this, we ex- 
plored the velocity and temperature proflle of the reversed stagnation-point flow. 
In present studies, investigations on the behaviors of dimensionless velocity in the 
reversed stagnation-point flow reveal that: 

1. Compared to the previous research, it is not quite appropriate to say Proud- 
man and Johnson are wrong because of neglecting the viscous term in their 
analytic result for region sufl&cient far from the wall. Also, their inviscid 
result is impressing; because one can expect the flow pattern (see Figure 



3.1) from the inviscid fleld. 



2. In the region near the reversed stagnation point the numerical simulation 
agrees remarkably well for smaller values of r with the known similarity 
solution, thus conflrming the predictions of the viscous Proudman- Johnson 
solution. Their idealization that the viscous forces are of the same order as 
the inertial forces is acceptable near the stagnation point. 

3. Separation will occur near the wall as 77 ^ and the region of reversed flow 
will move outward away from the wall. For large times r ^ oc, the reversed 
flow comes to rest. Viscous forces are dominant to decelerate the velocities 
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to zero and ultimately the region of reversed flow does not continue to grow 
but has finite dimensions. 

4. For the external flow outside the boundary layer, the hypothesis that the 
velocity v^x^rj^r) should pass over smoothly into that for inviscid Vq is 
not valid. The influence of backflow must be taken into account and a 
continuous change as Vo(r) decreases in magnitude for large r should be 
expected outside the boundary. 

On the other hand, investigations on the behaviors of dimensionless temperature 
in the nonisothermal reversed stagnation-point flow illustrate that: 

1. The solution of the thermal energy equation is also provided. The temper- 
ature distribution is monotonically increasing. The nondimensional tem- 
perature drops from its remote value to its wall value in a thin thermal 
boundary layer adjacent to the wall. It is surprising that, near the backflow 
region, there is no discernible temperature signature between the dividing 
streamlines. 

2. Larger Prandtl number results in thinner boundary layer and higher tem- 
perature gradient near the wall. When Prandtl number is small, the heat 
diffuses very quickly compared to the velocity field. This implies that for 
liquid metals the thickness of the thermal boundary layer is much bigger 
than that of the velocity boundary layer. 

3. The numerical simulation indicates that heat transfers mainly from the 
cooled wall to the heated external flow by conduction in the beginning. 
The heated flow passes though the wall, rises and cools down by conduction 
and convection of heat. Because of the motion of backflow, heated flow 
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sinks to the wall where it is prohibited from sinking further and becomes 
dense near the reversed stagnation point along the edges of the wall. 

With the establishment of this frame work, a similarity method applied to the 
two-dimensional unsteady reversed stagnation-point has induced new physically 
significant solutions, and application of the method to other case may be even 
more fruitful. Recommendations on the study of this type of fluid flow problem 
are given below: 

1. The similarity solution is valid only at the reversed stagnation point x — 0. 
In order to study the flow for non-zero values of x, we must revert the whole 
problem to the full boundary-equation. 

2. Three-dimensional simulation is much better than the two-dimensional case 
that we have been studying so far. However, more realistic simulation 
comes with high requirements in memory and CPU time so that the three- 
dimensional case is generally not simulated. A rapid development of com- 
puter hardware and software will further increase the opportunities for nu- 
merical simulation. 

3. More execution time would be sufficient in the simulation. Because of the 
time constraints, only a few cases of simulation are completed. More cases 
of simulations should be performed to obtain a more reliable data set of this 
type of fluid flow problem. 

4. In the result of numerical simulation, one may be observed that there 
are small vortices generated near the reversed stagnation point when the 
Reynolds is sufficient high. Some factors of affecting the probability of get- 
ting firm results of the investigations on the small vortices near the plate 
are thought to be: 
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(a) Sizes of the time steps; 

(b) Sizes of finite volume near the reversed stagnation point; 

(c) Magnitudes of the external flow velocity uq 

(d) Differences between the wall temperature and the ambient temper- 
ature Too 

5. The more important practical properties in engineering and technology ap- 
plication, like the velocity of wall is function of time r and the temperature 
of wall is function of time r, can be investigated and should be performed 
in the next phase of this study. 



APPENDIX 



MATLAB 

Listing 7.1: Finite-difference formulations for reversed stagnation-point flow 

1 clear all 

2 deta = 0.1; dtau = 0.05; Pr = l; 

3 IMAX = 100; NMAX = 6; 

4 beta=dtau / ( deta " 2) ; 

5 IM=IMAX/deta + l; NM=NMAX/dtau + l; 

6 h=zeros ( IM , NM) ; g=zeros ( IM , NM) ; s=zeros ( 1 , IM — 2) ; 

7 % 

8 % Setting the initial and boundary conditions 

9 eta(l ,1) = 0.0; 

10 for i= 2:IM 

11 eta(i,l) = eta(i — 1,1) + deta; 

12 h(i ,l)=l-erf (eta(i ,1) /2/ sqrt ( dtau/ 10) ) ; 

13 g(i,l)=erf(eta(i,l) /2/sqrt (dtau/lO/Pr) ) ; 

14 end 

15 

16 for n= 1 : NM 

17 h(l,n)=l; % no slip boundary condition 

18 g(lM,ii)=l; 
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19 end 

20 tic 

21 for n = l:NM— 1 



22 s(l ,l)=h(l ,ii)+0.5*beta*(h(3 ,ii)-2*h(2 ,ii)+h(l ,ii) ) ... 

23 +dtau*(2*h(l ,ii)-(h(l ,n) ) ^2) ... 

24 -dtau/2*(-h(3 ,ii)+4*h(2 ,ii)-3*h(l ,ii) ) * . . . 

25 0.5*(l-h(l ,ii))+0.5*beta*h(l ; 

26 for i = 2:IM-2 

27 s(l ,i)=h(i ,ii)+0.5*beta*(h(i + l,ii)-2*h(i , ii)+h ( i -1 ) 

28 +dtau*(2*h(i ,ii)-(h(i,ii))^2)... 

29 -dtau/2*(h(i + l,ii)-h(i-l,ii) ) * . . . 

30 (sum(l-h(l: i ) -0.5*(1 -h ( i , n) ) ) ; 

31 end 

32 s(l ,IM-2)=s(l ,IM-2)+0.5*beta*h(lM-2,ii) ; 

33 



34 % 

35 % Thomas algorithm for a tridiagonal system 

36 % a,b,c: diagonal , superdiagonal , 



37 % and subdiagonal elements 

38 a=(l+beta) *oiies ( 1 , IM — 2) ; 

39 b =— 0.5* bet a* ones ( 1 , IM — 2) ; 

40 c =— 0.5* bet a* ones ( 1 , IM — 2) ; 

41 x=ones ( 1 , IM — 2) ; 

42 d=ones ( 1 , IM — 2) ; 
d(l,l)=b(l,l)/a(l,l) ; 



y=ones(l ,IM-2) ; 
y(l,l)=s(l,l)/a(l,l); 

for p = l:(lM-3) 

den=a(l ,p + l)-c(l ,p + l)*d(l ,p) ; 
d(l ,p + l)=b(l ,p + l)/den; 

y(l,p + l) = (s(l,p + l)-c(l,p + l)*y(l,p))/den 

end 

x(l,IM-2)=y(l,IM-2); 
for p=IM-3:-l:l 

x(l,p)=y(l,p)-d(l,p)*x(l,p + l); 

end 

for p = 2:IM-l 

h(p,n+l)=x(l ,p-l); 

end 

s(l ,l)=g(l ,n)-dtau*(-g(3,n)+4*g(2,n)-3*g(l , 

*sum(l — h(l ,n) )+beta/Pr*g ( 1 ,n) ; 
for i = 2:IM-2 

s(l ,i)=g(i ,n) . . . 

-dtau/2*(g(i + l,n)-g(i-l,n)) 
(sum(l-h(l:i,n))-0.5*(l-h(i,n))) ; 

end 

s(l ,IM-2)=s(l ,IM-2)+beta/Pr*g(lM-2,n) ; 



70 % 

71 % Thomas algorithm for a tridiagonal system 

72 % a,b,c: diagonal , superdiagonal , 



73 % and subdiagonal elements 

74 a=(l+2*beta/Pr ) *ones ( 1 , IM-2) ; 

75 b=-beta/Pr*ones (1 , IM-2) ; 

76 c=— beta/ Pr* ones ( 1 , IM — 2) ; 

77 x=ones ( 1 , IM — 2) ; 

78 d=ones ( 1 , IM — 2) ; 

79 d(l,l)=b(l,l)/a(l,l) ; 

80 y=ones ( 1 , IM — 2) ; 
y(l,l)=s(l,l)/a(l,l) ; 

82 for p = l:(lM-3) 

83 den=a(l ,p + l)— c (1 ,p + l)*d(l ,p) ; 

84 d(l ,p + l)=b(l ,p + l)/den; 

85 y(l ,p + l) = (s(l ,p + l)-c(l ,p + l)*y(l ,p))/den; 

86 end 

87 x(l ,IM-2)=y(l ,IM-2) ; 

88 for p=IM-3:-l:l 

89 x(l ,p)=y(l ,p)-d(l ,p)*x(l ,p + l) ; 

90 end 

91 % — 

92 for p = 2:IM— 1 

93 g(p,n+l)=x(l ,p-l) ; 

94 end 



95 end 
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OpenFOAM 



Listing 7.2: myicoFoam solver 



1 / * 

2 Application 

3 myicoFoam 

4 

5 Description 

6 Transient solver for incompressible , laminar flow 

7 of Newtonian fluids and temperature profile 

8 \* 

9 

10 #include "fvCFD.H" 

11 

12 int main(int argc , char *argv[]) 

13 { 

14 ^include " setRootCase . H" 

15 ^include " createTime . H" 

16 ^include " createMesh . H" 

17 ^include " createFields . H" 

18 T^include " initCont inuityErrs . H" 

19 

20 Info« "\nStarting time loop\n" « endl ; 

21 



22 while ( runTime . loop ) 

23 { 
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24 Info«"Time = "«runTime . t imeName () « nl « endl ; 

25 

26 #iiiclude "readPISOControls .H" 

27 T^include " CourantNo . H" 

28 

29 f vVect orMatr ix UEqn 

30 ( 

31 fvm::ddt(U) 

32 + f vm : : div (phi , U) 

33 — f vm : : laplacian (nu , U) 

34 ) ; 

35 

36 solve (UEqn = — f vc : : grad (p) ) ; 

37 

38 // PISO loop 

39 

40 for (int corr = 0; corr<iiCorr ; corr++) 

{ 

42 volScalarField rUA = 1 . 0/ UEqn . A ( ) ; 

43 

44 U = rUA*UEqn.H() ; 

45 phi = ( fvc :: interpolate (U) & mesh.Sf()) 

46 + f vc : : ddtPhiCorr (rUA , U, phi); 

47 

48 adjustPhi (phi , U, p) ; 

49 
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50 for (int nonOrth = 0; nonOrth<=nNonOrthCorr ; 

nonOrth++) 

{ 

52 f vScalarMatrix pEqn 

53 ( 

54 f vm : : laplacian (rUA , p)= f vc : : div (phi ) 

55 ) ; 

56 

57 pEqn . set Reference (pRef Cell , pRef Value ) ; 

58 pEqn . solve ( ) ; 

59 

60 if (nonOrth = nNonOrthCorr ) 

{ 

62 phi — = pEqn.fluxO ; 

63 } 

64 } 
65 

66 T^include " cont inuityErrs . H" 

67 

68 U — = rUA*fvc : : grad(p) ; 

69 U. correctBoundaryCondit ions () ; 

70 } 
71 

72 // Temperature transport 

73 

74 f vScalarMatrix TEqn 



( 



f vm : : ddt (T) 

+ f vm : : div (phi , T) 

— f vm : : laplacian (DT , T) 



TEqn . solve ( ) ; 
runTime . write () ; 



Info« " Execut ionTime = " « runTime. 
elapsedCpuTime « " s" 
« " ClockTime = " « runTime. 

elapsedClockTime ( ) « " s" 
« nl « endl ; 



Info« "End\n" « endl; 
return 0; 
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Listing 7.3: Geometry Analysis 

/ \ 

1 FoamFile 

2 { 

3 version 2.0; 

4 format ascii ; 

5 class dictionary; 

6 object blockMeshDict; 

^ } 

8 convertToMeters 1; 

9 vertices 



10 ( 



11 


(0 





0) 


12 


(1 





0) 


13 


(1 


1 


0) 


14 


(0 


1 


0) 


15 


(0 





1) 


16 


(1 





1) 


17 


(1 


1 


1) 


18 


(0 


1 


1) 



19 ) ; 



20 blocks 

21 ( 

22 hex (0 1 2 3 4 5 6 7) (200 400 1) 

23 simpleGrading (1 5 1) 

24 ) ; 

25 edges ( ) ; 
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26 patches 

27 ( 

28 patch left_inlet ((2 6 5 1)) 

29 patch right.inlet ((0 4 7 3)) 

30 patch outlet ((3 7 6 2)) 

31 wall fixedWalls ((1 5 4 0)) 

32 empty front AndBack ((0 3 2 1) 

33 (4 5 6 7)) 

34 ) ; 

35 mergePat chPairs () ; 

V ) 
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Listing 7.4: Fluid Transport Properties 

/ \ 

2 FoamFile 

3 { 

4 version 2.0; 

5 format ascii ; 

6 class dictionary; 

7 location "constant"; 

8 object transportPropert ies ; 

^ } 

10 nu nu [ 2 -1 ] 2.94e-7; 

11 DT DT [ 2 -1 ] 1.68e-7; 

V J 
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Listing 7.5: Initial Pressure Profile 



FoamFile 



{ 



version 2.0; 

format ascii ; 

class volScalarField ; 

object p; 



10 



dimensions [02-20000]; 
internalField uniform 0; 



11 boundary-Field 



12 



13 



14 



15 



16 



17 



18 



19 



20 



{ 



} 



left_inlet {"type 
right_inlet {type 
outlet {"type 

value 
fixedWalls {"type 
f ront AndBack {type 



zeroGradient ; } 
zeroGradient ; } 
f ixedValue ; 
uniform 0;} 
zeroGradient ; } 
empty;} 
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Listing 7.6: Initial Velocities Profile 



FoamFile 

{ 

version 2.0; 

ascii ; 

volVectorField ; 
U; 



format 
class 
ob j ect 



} 



10 



dimensions [0 1 -1 0] ; 

internalField uniform (0 0); 



11 boundary-Field 



12 



13 



14 



15 



16 



17 



18 



19 



20 



21 



22 



{ 



left_inlet {"type 
value 

right_inlet {type 
value 

outlet {type 

fixedWalls {"type 



f ixedValue ; 
uniform (-1 0);} 
f ixedValue ; 
uniform (10 0);} 
zeroGradient ; } 
f ixedValue ; 



// 



value uniform (0 0);} 
frontAndBack {type empty;} 



// 
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Listing 7.7: Initial Temperature Profile 


1 


// 








2 


FoamFile 






3 


{ 








4 




version 


2.0; 




5 




format 


ascii ; 




6 




class 


volScalarField ; 


7 
8 


} 


ob j ect 


T; 




9 


dimensions 


[0 1 


0]; 


10 


internalField 


uniform 


373; 


11 


boundary-Field 






12 


{ 








13 




lef t _inlet 


{type 


f ixedValue ; 


14 






value 


uniform 373;} 


15 




right _inlet 


{type 


f ixedValue ; 


16 






value 


uniform 373;} 


17 




outlet 


{type 


zeroGradient ; } 


18 




f ixedWalls 


{type 


f ixedValue ; 


19 






value 


uniform 273;} 


20 




f ront AndBack {type 


empty;} 


21 


} 








22 


// 
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Listing 7.8: Simulation Control 



2 FoamFile 

3 { 

4 version 

5 format 

6 class 

7 location 

8 object 

^ } 

10 application 

11 startFrom 

12 startTime 

13 stopAt 

14 endTime 

15 delta! 

16 writeControl 

17 writ eint erval 

18 purgeWrite 

19 writeFormat 



2.0; 
ascii ; 
dictionary ; 
" system" ; 
controlDict 

myicoFoam ; 
StartTime ; 

0; 

endTime ; 

10; 

0.001; 
timeStep ; 
100; 

0; 

ascii ; 



20 writ ePrecision 6; 

21 wr it eCompression uncompressed; 

22 timeFormat general; 

23 t imePrecision 6; 

24 runTimeModif iable yes; 
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Listing 7.9: System Solver 



FoamFile 



{ 



version 

format 

class 



} 



10 



11 



12 



13 



14 



15 



16 



17 



18 



19 



20 



21 



22 



23 



24 



25 



2.0; 
ascii ; 
dictionary ; 



location "system" ; 
object fvSchemes 



ddtSchemes 

{ 

default 

} 

gradSchemes 

{ 

default 

grad(p) 

} 

divSchemes 
{ 

default 

div(phi ,U) 
div(phi ,T) 

} 

laplacianSchemes 



Euler ; 



Gauss linear; 
Gauss linear; 



none ; 

Gauss linear; 
Gauss upwind; 
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26 { 



27 



default none ; 

28 laplacian (nu , U) Gauss linear corrected; 

29 laplacian (( 1 I A(U) ) ,p) Gauss linear corrected; 

30 laplacian (DT , T) Gauss linear corrected; 

31 } 

32 int erpo 1 at ionSc hemes 

33 { 

34 default linear; 

35 int erpolat e ( HbyA ) linear; 

36 } 

37 snGradSchemes 

38 { 

39 default corrected; 

40 } 

41 fluxRequired 

42 { 

43 default no ; 

44 p ; 

45 } 
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Listing 7.10: Preconditioner and Tolerance 



FoamFile 



{ 



version 
format 
class 
location 
ob j ect 



2.0; 
ascii ; 
dictionary ; 
" system" ; 
f vSolut ion ; 



} 



10 



11 



12 



13 



14 



15 



16 



17 



18 



19 



20 



21 



22 



23 



24 



solvers 

{ 

p 

{ 



solver 



PCG 



preconditioner DIG; 
tolerance le— 06; 



relTol 



0; 



} 
T 

{ 



solver 
preconditioner 
tolerance 
relTol 



BICCG ; 
DILU ; 

le-7; 
0; 



25 
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26 U 

{ 

28 solver PBiCG; 

29 precondit ioner DILU; 

30 tolerance le— 05; 

31 relTol 0; 

32 } 



33 } 

34 PISO 



35 { 

36 nCorrectors 2; 

37 nNonOrthogonalCorrect ors 0; 

38 pRefCell 0; 

39 pRefValue 0; 

40 } 

V / 
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SPECIFICATIONS OF THE SIMULATION COMPUTER 



Computer Model 


Lenovo Thinkstation Workstation D20 


CPU 


Intel ®Xeon ©CPU X5690 @3.47 GHz 


RAM 


24.0 GB 


Operation System 


Ubuntu Linux 10.04 
Windows 7 


Software 


OpenFOAM 
MATLAB 
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